Elasticity, fluctuations and vortex pinning in ferromagnetic superconductors: 

A "columnar elastic glass" 

A. M. Ettouhami 1 , Karl Saunders 2 , L. Radzihovsky 3 and John Toner 4 
1 Department of Physics, University of Florida, P.O. Box 118440, Gainesville, FL 32611-8440 
2 Department of Physics, California Polytechnic State University, San Luis Obispo, CA934 07 
3 Department of Physics, University of Colorado, 390UCB, Boulder, CO 80309-0390 
4 Department of Physics, Materials Science Institute and Institute of Theoretical Science, 
University of Oregon, Eugene, OR 97403 
(February 2, 2008) 

We study the elasticity, fluctuations and pinning of a putative spontaneous vortex solid in ferro- 
magnetic superconductors. Using a rigorous thermodynamic argument, we show that in the idealized 
case of vanishing crystalline pinning anisotropy the long-wavelength tilt modulus of such a vortex 
solid vanishes identically, as guaranteed by the underlying rotational invariance. The vanishing of 
the tilt modulus means that, to lowest order, the associated tension elasticity is replaced by the 
softer, curvature elasticity. The effect of this is to make the spontaneous vortex solid qualitatively 
more susceptible to the disordering effects of thermal fluctuations and random pinning. We study 
these effects, taking into account the nonlinear elasticity, that, in three dimensions, is important at 
sufficiently long length scales, and showing that a "columnar elastic glass" phase of vortices results. 
This phase is controlled by a previously unstudied zero-temperature fixed point and it is charac- 
terized by elastic moduli that have universal strong wave- vector dependence out to arbitrarily long 
length scales, leading to non-Hookean elasticity. We argue that, although translationally disordered 
for weak disorder, the columnar elastic glass is stable against the proliferation of dislocations and 
is therefore a topologically ordered elastic glass. As a result, the phenomenology of the sponta- 
neous vortex state of isotropic magnetic superconductors differs qualitatively from a conventional, 
external-field-induced mixed state. For example, for weak external fields H, the magnetic induction 
scales universally like B{H) ~ B(0) + cH a , with a « 0.72. 



I. INTRODUCTION 

It has been more than twenty years since supercon- 
ducting materials containing a periodic arrangement of 
magnetic rare earth ions were discovered. 1 These in- 
cluded the rare-earth compounds i?Rh4B4, R^MogSs and 
RajMogSeg, with x = 1.0 or 1.2 and R a rare earth mate- 
rial (such as Lu, Y, Tm, Er, Ho and Dy). As temperature 
was varied, these materials displayed either a supercon- 
ducting or a magnetically ordered phase, but no phase 
where both superconductivity and magnetism coexisted 
simultaneously. More recent efforts, driven by the need to 
understand high-T c superconductivity and other strongly 
correlated materials, have led to experimental discov- 
eries of systems exhibiting coexistence of superconduc- 
tivity and magnetic order, such as the rare-earth nickel 
borocarbides 2 i?Ni2B2C with relatively high supercon- 
ducting transition temperatures (T c ~ 16.5 K for the Lu 
and 15K for the Y compound). These materials exhibit 
a rich phase diagram that includes superconductivity, 
antiferromagnetism, ferromagnetism and spiral magnetic 
order. 3,4 While most of these materials are antiferromag- 
nets, there is now ample experimental evidence 5 that be- 
low both the Curie (~ 2.3 K) and the superconducting 
transition temperatures (~ 10.5 K), superconductivity 
and ferromagnetism competitively coexist in ErN^I^C 



compounds. Other possible examples of ferromagnetic 
superconductors (FS) are the recently discovered high 
temperature superconductor S^YRui-^CuzOe and the 
putative p-wave triplet strontium ruthenate supercon- 
ductor, Sr2Ru04, which spontaneously breaks time re- 
versal symmetry, as well as the recently discovered com- 
pounds RuSr 2 GdCu 2 8 , 6 UGe 2 , 7 ZrZn 2 , 8 and URhGe. 9 

While there are many interesting unanswered micro- 
scopic questions regarding the nature of such ferromag- 
netic superconductors, much of their low-energy phe- 
nomenology is dictated by general symmetry principles. 
These constrain the form of the Ginzburg-Landau the- 
ory involving the local magnetization M and the super- 
conducting order parameter 'J', that we expect to de- 
scribe the low-energy equilibrium thermodynamics. As 
we discuss in Sec. II, the resulting theory predicts a 
rich phenomenology, 10-16 that among other phases ad- 
mits a very interesting spontaneous vortex (SV) state 
driven by the spontaneous magnetization, rather than 
by an external magnetic field H. The argument made 
in the context of borocarbides 15,16 (which we believe ap- 
plies more generally) implies that ferromagnetic super- 
conductors are expected to exhibit such a spontaneous 
(H = 0) vortex state. Recent small angle neutron scat- 
tering experiments 17 on ErNi2 11 B2C provide some evi- 
dence for the existence of a SV state. The combined 
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effects of soft elasticity, random pinning, and thermal 
fluctuations lead to a unique phenomenology for the SV 
solid; this unique phenomenology is the subject of this 
paper. 

We will show, in particular, that for H = 0, and van- 
ishing crystal anisotropy pinning fields 18 the elastic prop- 
erties of the resulting SV solid differ dramatically and 
qualitatively from those of a conventional Abrikosov lat- 
tice. The key underlying difference is the vanishing of 
the tilt modulus, which is guaranteed by the underly- 
ing rotational invariance. 18 This invariance is of course 
broken by the magnetization M, but because this sym- 
metry breaking is spontaneous, the tilt modulus remains 
zero. This is reflected by the invariance of the energy un- 
der a simultaneous global rotation of the magnetization 
and vortex lattice induced by it. Formally, the vanish- 
ing of the tilt modulus corresponds to the vanishing of 
the "mass" of the local field which is a Goldstone mode 
associated with the rotational symmetry that is sponta- 
neously broken by the vortex lattice and magnetization. 
This contrasts strongly with a conventional vortex solid, 
where the rotational symmetry is explicitly broken by the 
applied field H and crystalline anisotropy. Of course, in 
a real sample there will always be some degree of crys- 
talline anisotropy. We will show that our results for the 
columnar elastic glass phase of the SV solid are valid out 
to a long length scale that depends on the strength of 
crystalline anisotropy pinning fields. This length scale 
diverges with vanishing anisotropy so that our results 
best describe materials with weak crystalline anisotropy. 

All of the conclusions that we draw about other dis- 
tinctive properties of the SV solid for H = are a direct 
consequence of this important observation. In particular, 
we find that this "softness" (i.e., vanishing tilt modu- 
lus) of the SV lattice dramatically enhances the effects of 
quenched disorder and thermal fluctuations. As in con- 
ventional vortex lattices, 19 any amount of disorder A, 
however weak, is sufficient to destroy translational order 
in SV lattices. Here we find that the Larkin lengths Rjr 
and R z c , beyond which translational order is destroyed by 
random forces and torques, are highly anisotropic both in 
their magnitude and in the way they scale with pinning 
strength A. We predict 

i?^ocl/A 2 / 3 , (1.1a) 
^ocl/A 1 / 3 , (1.1b) 

in d = 3 dimensions, a prediction that can in principle 
be checked in neutron scattering and transport measure- 
ments. However, as we will show, for sufficiently weak 
pinning (such that dislocation loops remain bound and 
the elastic description remains valid) the disorder also 
qualitatively alters the SV solid's long-scale elastic prop- 
erties, leading to "anomalous elasticity": the universal 
scaling of elastic moduli with wavevector q out to ar- 
bitrarily long length scales (small q), with some elastic 
moduli vanishing and others diverging as the wavevector 
q — > 0. Other related long length scale elastic prop- 



erties include a universal negative Poisson ratio and a 
non-Hookcan elasticity. This anomalous elasticity is com- 
pletely distinct from the well-known wavevector depen- 
dence of the tilt and bulk moduli at short length scales 
(less than the London penetration length A, or, equiv- 
alently, at wavevectors q such that qX 3> 1) in con- 
ventional vortex lattices. 20 This long-wavelength anoma- 
lous behsvior is characteristic of a new kind of topo- 
logically ordered "columnar elastic glass" phase of vor- 
tices, which can be shown 21 to be stable, for weak disor- 
der, against the proliferation of dislocations. This spon- 
taneous vortex solid is another example of a class of 
"soft" elastic systems 22 , where long-scale elasticity of the 
strongly-interacting Goldstone modes (associated with 
the spontaneously-ordered phase) is characterized by a 
nontrivial low-temperature fixed point leading to univer- 
sal anomalous elasticity. 

The best way to experimentally probe the conse- 
quences of spontaneously broken symmetry is to break it 
directly with a weak external field. Here we predict that, 
as a consequence of the anomalous elasticity, the increase 
in the magnetic induction SB(H) = B(H) — B(0) over 
the spontaneous induction B(0) as a function of weak 
applied field H obeys a universal "non-Hookean" scaling 
law: 

SB(H) oc H a , (1.2) 

with the universal exponent a = 0.72 ± 0.04 (see Fig. 1), 
a prediction that should be experimentally testable. This 
non-Hookean scaling is in contrast to conventional vor- 
tex lattices in which SB(H) scales approximately linearly 
with H - H cl for H > H cl . 28 

The rest of this article is organized as follows. In Sec. 
II, we review the Ginzburg-Landau theory of magnetic 
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FIG. 1. The nonlinear and universal power-law B(H) 
scaling, which at weakest fields H < H cr , and strongest 
fields H > Hnl is cutoflt by the crystal symmetry breaking 
anisotropy and £ NL , respectively. 



2 



superconductors, considering the idealized case of van- 
ishing crystalline pinning anisotropy, 18 and use it to de- 
rive the effective interaction between vortex lines in such 
materials. Using this result, we show rigorously, using 
a thermodynamic argument, that as a consequence of 
the underlying rotational invariance the tilt modulus at 
long wavelengths ca{q z = 0) vanishes in the spontaneous 
vortex state, thereby showing that the low-energy, long- 
wavelength tilt modes of the spontaneous vortex lattice 
are controlled by curvature elasticity of the form k(i9 2 u) 2 
(rather than by a conventional tension type (<9 z ii) 2 , with 
u(r) the transverse displacement field of the vortex lat- 
tice), and derive the curvature modulus n. In Sec. Ill, 
we examine the consequences of the softness of the tilt 
modes for the physics of the fluctuations and pinning of 
the spontaneous vortex lattice. We show that while both 
are strongly enhanced quantitatively, in clean samples we 
expect the spontaneous vortex lattice phenomenology to 
remain unmodified. This is because the thermal mean- 
squared fluctuations in the SV solid, in three dimensions, 
are dominated by contributions from short length scales 
(at which the spontaneous lattice elasticity should not 
differ qualitatively from that of a conventional lattice.) 
In particular, we expect only quantitative changes in the 
melting curve of a spontaneous vortex lattice in a clean 
sample. In contrast, we find that the effects of pinning 
by random impurities differ qualitatively from that of a 
conventional vortex lattice. Specifically, we show that 
for d < d v uc = 9/2, (in contrast to d v uc — 4 in an ordi- 
nary vortex lattice) random pinning disorders the spon- 
taneous vortex lattice on length scales larger than the 
highly anisotropic Larkin domains given by Eqs. (1.1a)- 
(1.1b) (sec also Eqs. (4.14a)-(4.14b) below). In Sec. IV 
we begin by ignoring the effects of nonlinear elasticity 
and tilt disorder, and show that on length scales longer 
than R z c ' , tne power-law roughness of the vortex lattice 
crosses over to the logarithmic growth of mean-squared 
phonon distortions, as in conventional vortex lattices. 
Next, we include tilt disorder (that is always present) 
and find that for d < = 7/2 it leads to a far stronger 
power-law roughness of the vortex lattice, thereby domi- 
nating over the logarithmic roughness caused by the vor- 
tex positional pinning. In Sec. V, we show that as a 
consequence of such strong power-law distortions, the 
nonlinear elasticity becomes qualitatively important, as 
indicated by the breakdown of perturbation theory in 
the presence of an (in principle) arbitrarily weak pin- 
ning. These findings are consistent with similar results 
for other soft solids in the presence of disorder, most 
notably from studies of liquid crystals confined in a ran- 
dom matrix 22 ' 29 . In Sec. VI we use the momentum-shell 
renormalization group method to determine the conse- 
quences of these effects for the long-scale properties of 
the randomly-pinned spontaneous vortex lattice. We 
show that for scales larger than an anisotropic nonlin- 
ear crossover length scale, the elasticity becomes 
anomalous and scales universally with wavevector. It is 
controlled by a nontrivial zero-temperature fixed point 



that is perturbative in e = 7/2 — d. In Sec. VII we 
analyze the topological stability of the SV solid and ar- 
gue that the resulting phase in the presence of disorder 
is a "columnar elastic glass" 30 ' 21 (CEG) that is stable 
against the proliferation of dislocations. The term "elas- 
tic glass" describes a phase that, while translationally 
disordered, retains topological order. Other examples of 
elastic glasses have been predicted for smectic 29 ' 31 ' 32 and 
columnar (discotic) 21,33 liquid crystals. Associated with 
the loss of translational order is a broadening of the Bragg 
scattering peaks, so that they now exhibit power-law de- 
cay rather than delta-function divergences. The absence 
of even quasi-sharp Bragg peaks in these systems is the 
reason we refer to them as "elastic glasses" rather than 
Bragg glasses (of which the ordinary (field-induced) dis- 
ordered vortex lattice is an example 34 ). An elastic glass 
is distinguished from the liquid phase by being free of 
unbound dislocation loops, which proliferate in the liq- 
uid. The glass transition is identified as an "unbinding" 
of the dislocation loops. This transition is then, quali- 
tatively, very similar to the melting of the flux lattice in 
the absence of disorder, which can also be thought of as 
an unbinding of dislocation loops. The only difference is 
that in the glass problem, the flux lattice is translation- 
ally disordered both above and below the transition. The 
absence of defects below the transition, however, means 
that the low temperature solid phase still has a finite 
shear modulus, leading to glassy behavior. The exper- 
imental consequences of this behavior are discussed in 
Sec. VIII, and Sec. IX contains a discussion of possible 
future theoretical directions. 



II. GINZBURG-LANDAU PHENOMENOLOGY 
AND HARMONIC ELASTIC THEORY OF A 
SPONTANEOUS VORTEX LATTICE 

Irrespective of the microscopic mechanism that is ul- 
timately responsible for the stability of the spontaneous 
ferromagnetic superconducting state, general symmetry 
principles dictate that the long length scale phenomenol- 
ogy is described by an effective Landau-Ginzburg free 
energy functional 

F^/, r {I a( T)H 2 + J|^ + |l|(V-^a),M 2 + 
+ 7T + 7: a (T) |M| 2 + - /? |M| 4 + - 7 1 VM| 2 + 

07T 2 4 2 

b-M — b-H}, (2.1) 

47T J 

where ip(r) is the local superconducting order parame- 
ter, M(r) is the magnetization, <p — hc/2e is the quan- 
tum of flux, a is the vector potential, b = V x a is the 
total flux density inside the superconductor, H is the 
external magnetic field, and we have (for now) ignored 
crystalline anisotropy. 18 The constants a, b, a, [3 and 7 
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are experimentally accessible phenomenological param- 
eters, with a(T) and a = a(T) — Att changing sign at 
the superconducting and ferromagnetic transition tem- 
peratures, respectively. 35 The magnetic field h associated 
with the superconducting currents satisfies Ampere's law 
V x h = 4-7rj/c (c is the speed of light in vacuum) and is 
related to b by: 



b(r) = h(r) + 47rM(r) + H. 



(2.2) 



The term |VM| 2 in Eq. (2.1) stands for E*( VM 2 - 
From Eq. (2.2) we observe that the flux density b is de- 
termined by magnetization of the local moments M, the 
external magnetic field H and the screening magnetic 
field h generated by the diamagnctic orbital currents j. 

The minimization of F with respect to tp, a and M 
leads to the following three coupled equations: 



2m 
1 



ij- a) V = ~ aip - \ b\il>\ 2 i) 



— ,— \ ieh 
— V x (V x a) = 

4tt v ; 2mc 



2^2 
ipVtp* - ip*Vtp 



+ 



2e 2 1 

'•' 2 a + V x M + —V x H , 



47T 

b(r) = aM + /3|M| 2 M - 7 V 2 M . 



(2.3a) 



(2.3b) 
(2.3c) 



It has been shown many years ago by Kuper et al. 12 and 
by Grccnside et al. 14 that, depending on the actual val- 
ues of model parameters, the above mean-field equations 
admit a variety of equilibrium phases, with transitions 
between them controlled by the temperature, which en- 
ters model parameters, most notably a(T) and a(T). In 
addition to the paramagnetic normal phase, M = ip = 
appearing for a(T) > 0, a(T) > 0, a pure superconduct- 
ing state, M — 0, ip ^ 0, a ferromagnetic normal state 
(M 7^ 0, i\) = 0) and a stable spiral phase (where super- 
conductivity coexists with spiraling magnetization) can 
appear. Furthermore, the Ginzburg-Landau theory also 
predicts a thermodynamically stable spontaneous vortex 
state for a range of physically realistic parameters that 
give a large Abrikosov ratio required for a robust mixed 
state (A/£ <~ O(10), £ being the superconducting coher- 
ence length and A the London penetration depth), and 
a large exchange 7 necessary to suppress the competing 
spiral phase. In this mixed state vortices are induced 
by the local magnetic moments that are spontaneously 
ordered into a uniform ferromagnetic state. This is in 
contrast to conventional vortex lattices that are induced 
by an external magnetic field. 

It is not our intention in this article to study the rich 
phase diagram that follows from Eqs. (2.3a)-(2.3c). For 
this, we refer the interested reader to previous work 12 ' 14 , 
where this has been done in great detail. Instead our fo- 
cus here is on the spontaneous vortex state. We use the 
above Ginzburg-Landau theory to derive the vortex line 
interaction and elasticity of the resulting spontaneous 
vortex solid phase and study it in the presence of thermal 
fluctuations and random pinning. 



A. Interaction potential between vortices in the SV 

state 

The description of the vortex state can, as usual, be 
derived from the Ginzburg-Landau theory by reexpress- 
ing the fields ip and a in terms of the vortex positions, 
thereby reexpressing energy in terms of the vortex con- 
formational degrees of freedom. To this end, taking the 
curl of equation (2.2), using Ampere's law V x h = ^ZLj, 
and b = V x a, we obtain 

-3- V x (V x a) = -j(r) + V x M(r) (2.4) 

for zero external magnetic field. 

Comparing this expression with Eq. (2.3b) gives the 
standard expression for the superconducting current: 28 



ieh 
2m~ 



2e 2 



mc 



IVfa, (2.5) 



which, in the London approximation ip = ipoe ll) ( r \ (with 
ipo a constant except inside the small vortex core of size 
£, which is valid in the large A/£ ^> 1 limit), reduces to 28 



J<"~jef«'>-s v «- 



(2.6) 



This relates superconducting currents to the vector po- 
tential and the phase of the superconducting order pa- 
rameter, with the London penetration depth given by 

A = (mc 2 /47re 2 i/>g) 2 . Taking the curl of this current and 
using Ampere's law, together with Eq. (2.2) and 



V x V0 = 2ir^8 2 {r-v v {z)) 



(2.7) 



dz 

V 

where the sum runs over vortices, with 

R u {z) = (v v {z),z) (2.8) 

parametrizing the trajectory of the ^th vortex line as it 
traverses the superconducting sample (we choose the av- 
erage magnetization and vortex direction to be along the 
z axis), we find: 



b(r) - A 2 V 2 b(r) = faJ2 5 ^ r ~ r "( z » 

V 

- 47rA 2 V 2 M(r) . 



dz 



(2.9) 



This London equation differs from its counterpart in 
ordinary superconductors by the presence of the local 
magnetization term on its right hand side. It needs to be 
solved simultaneously with the constitutive relation, Eq. 
(2.3c), relating flux-density b and magnetization M. 

Using Eq. (2.3c) and the H = version of Eq. (2.2), 
we can write the following equation for the screening field 
h(r): 
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h(r) =5M + /3\M\ 2 M - 7 V 2 M 



(2.10) 



with a = a — Air. For small h(r) and a < this 
equation has the spatially uniform solution M(r) = 

h\J ^ = Mo. Linearizing about this solution by writing 

M(r) = Mo + £M(r) and expanding to linear order in 
<5M(r) gives, for q ^ 0, 



M(q) = <JM(q) = X (q)h(q), 



(2.11) 



where the second equality follows from the fact that Mo, 
being spatially uniform, contributes nothing to M(q) for 
q ^ , and 



x(q) = 



2\a{T)\+ lq 2 



(2.12) 



Combining this equation with Eq. (2.9) gives, in Fourier 
space: 



(1 + AV)b(q) = 0o J dze ~ 



iq-Rjy (z) 



dll v 

dz 



+ 4ttA Vx(q)Mq) . 



(2.13) 



Now, using Eq. (2.11) together with the H = version 
of Eq. (2.2) gives: 



b(q) = (1 + 47T X (q)) h(q) : 



(2.14a) 



which allows us to solve for the magnetic field h(q) and 
flux density b(q): 



-iq-R„(z) 



+ AV+47T X (q) ~dT' 



h(q)=0o^ dz T 

v J 

^ f (l + 4 7 r X (q))e-^ R -( z ) dH v 
b(q) = cfro ^Jdz 1 + Aag3 + 4wx(q) 



(2.15a) 
(2.15b) 



It is easy to verify that, for one flux line, / dx b z (x, z) = 
<fro, giving the magnetic flux quantization for a single vor- 
tex. On the other hand, the magnetic field h(r) due to 
the screening currents around the vortex cores satisfies 
J dxh z (x,z) ~ (fro/(l + 47rx(0)) 7^ (fro and is therefore 
not quantized. 

To proceed further we use the London approximation 
for the kinetic energy density 



h .,_ 2m . 
2^ I(V - fc" a > 



2ttA 2 2 
—2-3 , 



(2.16) 



to reduce the H = Ginzburg-Landau free energy func- 
tional (2.1) to the London expression, omitting a con- 
stant associated with the condensation energy: 



f , r h ,, 2in . ,,, b 2 

+ i a(T) |M| 2 + 1 /? |M| 4 + 1 7 | VM| 2 - b • m}, 



i a (T)|M| 2 + ^|M| 4 + i 7 |VM| 2 



b-M|. 
(2.17) 



Using the fact that b = h + 47rM, the London free energy 
reduces to a sum F L — F s + F m , where (we henceforth 

use the shorthand notation f q = / ) 



Sir J 



/i 2 (r) + A 2 (Vx h) 5 



2„2 



1 + A 2 g 
8tt 



|h(q) 



(2.18) 



is the free energy associated with superconducting cur- 
rents, and 



d A r 



i<5|M| 2 + i/3|M| 4 + i 7 |VM| 2 



L 



2X(q)|h(q)| 2 



(2.19) 



is the ferromagnetic part. We thus obtain: 
1 f rl + \ 2 q 2 



F L 



Air 



+ X(q)]|h(q)| 2 . (2.20) 



Now, using the expression (2.15a) for h(q) we finally ob- 
tain the London free energy of an arbitrary conformation 
of interacting flux lines: 

(2.21) 

with the interaction potential: 



V(R^z)-R„(z')) 



-»q-(R„(z)-R„(z')) 

4tt 7 q 1 + A V + 4^ X (q) ' 



(2.22) 



In the absence of the ferromagnetic component (x(l) — 
0), V(q) reduces to the usual London interaction between 
vortices in a conventional isotropic superconductor. 

Using our earlier result (2.12) for x(q) implies that the 
interaction potential V is given in Fourier space by: 



f(q) 



2\&\+ ig 2 



Att (1 + AV)(2|a| + iq 2 ) + An ' 



(2.23) 



Expressions (2.12) and (2.23) are the expressions we shall 
use below to derive the elastic moduli of the SV lattice. 
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B. Tilt and curvature moduli, and harmonic elastic 
theory of a spontaneous vortex lattice 

Having derived the interaction potential between flux 
lines, we are now in a position to calculate the tilt modu- 
lus of the spontaneous vortex lattice and to show that in 
the absence of an external field, it indeed vanishes identi- 
cally as dictated by rotational symmetry. 18 To this end, 
we shall write the vortex trajectories as 



t v {z) =X„ + u(X„, z), 



(2.24) 



with Xj, the equilibrium lattice position of the vth flux 
line and u(X„,z) the two component displacement at 
height z relative to X„. Expanding the London free en- 
ergy Fl of the spontaneous vortex system, Eq. (2.21), up 
to quadratic order in the displacements u(X y , z) leads to 
Fl = Fq + F e i , where Fo is the energy of an undistorted 
vortex lattice: 

F ° = \ E / dz I dz ' v ( x » - x «" z - z ') > 



W]Tv(Q). 



(2.25) 



Q 



and F e i is the elastic energy of the lattice, which in 
Fourier space is given by: 



1 



Ua(q)*a/3(q)w/?(-q) , 



(2.26) 



where Q is the volume of the system, n — B/<j>o is the 
average density of flux lines (B being the magnetic in- 
duction), the Q's are the reciprocal lattice vectors, and 
we have used Einstein's implicit summation convention 
over the repeated indices a and (3, which run over the 
two directions (a; and y) perpendicular to z, the mean 
direction of the flux lines. The elastic matrix $ Q/ g(q) of 
the FLL derived from using the expansion (2.24) in the 
clastic Hamiltonian (2.21) has the usual form: 20 

$a/3(q) = n 2 { [(<?« - q a )(Qf3 - qp) 
Q 

+ q%p]V(Q-q)-Q a QpV(Q)} (2.27) 

with V(q) the interaction potential of Eq. (2.23). 

As discussed in the Introduction, in a SV state with no 
external field to explicitly break rotational symmetry, 18 
the underlying rotational invariance guarantees a vanish- 
ing of the vortex- line tilt modulus C44. We now prove this 
explicitly for a spontaneous vortex lattice, via a rigorous 
thermodynamic argument. 

Given that the elastic matrix of a hexagonal vortex lat- 
tice (which is the lattice type we expect for an isotropic 
interaction potential) is of the general form 20 (in the fol- 
lowing, and throughout the rest of the paper, stands 
for the projection of q onto the plane that is orthogonal 



to z, and greek indices in q a , qp, . . . run over components 
of qj_ only): 

*a/?(q) = (cii - c 66 )q a qp + 5 a p(cuq 2 z + c 6e q 2 ± ) , (2.28) 

we see that the tilt modulus C44(q) can be extracted from 
<fr a p(q z z), which, according to (2.28), should be simply 



*a/3(g z z) = C44(q z z)q 2 z S a p. 



(2.29) 



Setting q = q z z in our general expression (2.27) for <& a p 
gives 



<t> a p(q) =n 2 J2 QcQp \V{Q - q z z) - V(Q) 
Q 

+ q 2 z S a pV(Q - q z z) 



(2.30) 



Using the fact that V(q) depends only on the squared 
magnitude q 2 of q , we can rewrite this as: 



<f a p(q z z) = n 2 Y J Q a Qp\V{Q 2 + ql)-V(Q 2 ) 
Q 

->2 1 „2\„2 



+ V(Q 2 +q 2 z )qj6 a p 



(2.31) 



where V is the interaction potential such that V(q) = 
V(q 2 ). Now, using the fact that for a hexagonal lattice 

E QaQpfiQ 2 ) = IJ2 Q 2 f(Q 2 ) s ^ , (2.32) 

Q Q 

we can rewrite (2.31) as 

® a p(q z z) = s a pn 2 ]T % \v(Q 2 + q 2 z )-v(Q 2 ) 



Q 

2 1 „2\„2 



+ V(Q z +q z z )q 



(2.33) 



which is exactly of the form required by (2.29), enabling 
us to identify 



C4i(q z ) = n 2 ^2 



2^Q 2 [V(Q 2 + q 2 z )~v(Q 2 ) 



+ V(Q 2 + q 2 z ) 



(2.34) 



This expression of c 4 4(q z ) has the following limit as 
q z -> 0: 

c 44 (0) = n 2 [ ^'(Q 2 ) + V(Q 2 ) ] , (2.35) 
Q 



which can be rewritten as: 20 

C44 (0Hi?A[in 2 ]>>(Q)]. 



(2.36) 



To see this last step more explicitly we first note that 
keeping the lattice structure hexagonal, increasing the 
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magnetic field B must decrease the lattice constant a 
(since the flux carried by each vortex line is a fixed 
quantum of flux (j) ). Indeed, since for hexagonal lat- 
tice B — 4>q I '(a 2 y/3) , we see that Q 2 oc l/a 2 (B) oc B, for 
all reciprocal lattice vectors Q. Equivalently, 

Q 2 (B) = Q 2 (B )B/B , (2.37) 

for all reciprocal lattice vectors Q. It then follows that 



0V(Q 2 ) = y, m dC? 



dB 



= V'(Q 2 ) 



dB 



Q 2 (B ) 

B n 



(2.38) 



where denotes derivatives with respect to B keep- 

ing the lattice structure fixed, and we have used (2.37) in 
the last two equalities. Now using this result (2.38) and 
n = B/4>o , we have 

Q Q 



(2.39) 



as claimed earlier. 

As can be seen from Eq. (2.25) above, the quantity in 
square brackets in this last equation is nothing but the 
free energy density (Fo/Q) of a lattice of straight vortex 
lines, which leads to: 



CAA{q z =V) = B—(F /n). 



(2.40) 



In the presence of an external magnetic field H, thermo- 
dynamics dictates that 



d(F /n) = H_ 
dB 4tt 



and hence we find 



ca{q z = 0) 



BH 



(2.41) 



(2.42) 



which shows, that for H = the long-wavelength tilt 
modulus vanishes identically 



c\a{1z = 0) = 0, 



(2.43) 



as required by rotational invariance. 

Having established the vanishing of the tilt modulus 
at q z — 0, let us now find the limiting form of C44(q z ) 
for small q z . Using a Taylor expansion of the right hand 
side of Eq. (2.34) near q z — leads to the following 
expression 



C44{qz )^n 2 q 2 z ^[^V"(Q 2 ) + V'(Q 2 
Q 



(2.44) 



Neglecting the periodicity of the flux lattice and retain- 
ing only the Q = term (which is usually dominant at 
high density) in the sum gives: 



C4i(q z )- 



B 2 



B 2 /7r 7 -5 2 A 2 \ 

With K ~ -— — — -T7T 

4tt V[27t+ |<5|1 2 / 



(2.45a) 
(2.45b) 



Hence, the long- wavelength bending energy of the spon- 
taneous vortex lattice is characterized by a curvature (as 
opposed to tension) elastic energy of the form k(<9 2 u) 2 . 
Note, however, that expression (2.45a) and expression 
(2.45b) for the curvature modulus k are valid only at 
long wavelengths, i.e. for q z — ► 0. At short wavelengths, 
we can restrict the reciprocal lattice sum in Eq. (2.34) 
to the Q = term (the so-called continuum approxi- 
mation), upon which the tilt modulus c^q) takes the 
following expression: 



c 44 (q) - J 



1 



+ X 2 q 2 + 4^ X (q) 



(2.46) 



where we used the fact that n — B/fo. In the absence 
of a magnetic component, x(q) = and the above ex- 
pression reduces to the usual result for the tilt modulus 
in conventional type-II superconductors. 

We can now combine this result for the tilt modulus 
with the standard analysis 20 usually applied to obtain 
the compression and shear moduli and thus write down 
the harmonic elastic energy for the spontaneous vortex 
lattice. From the general expression of the elastic matrix, 
Eq. (2.28), the compression and shear moduli are given 
by cu(q) = $ xx (qx)/q 2 and c e6 (q) = $xx(qy)/q 2 respec- 
tively. Hence, using Eq. (2.27) for $(q), we obtain: 



c n(<?) = ^ E I®* - ^)MQ - ?x)-Q 2 T/(Q)], 



q Q/0 



(2.47a) 



(2.47b) 



The compression modulus is readily evaluated in the con- 
tinuum (Q = 0) limit, with the result: 



cu(q) 



B 2 



1 



4tt 1 + A 2 9 2 + 4ttx(<7) ' 



(2.48) 



(where we have again used the fact that n — B/(j) ). On 
the other hand, in expression (2.47b) for the shear modu- 
lus, expanding V(Q — qy) to second order in q we obtain 
the following expression: 
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2 h dQ y 



(2.49) 



a result which should be valid (c@g being in general 
weakly dispersive) over the entire Brillouin zone of the 
SV lattice. These results are identical to the standard ex- 
pressions for the ordinary (field-induced) vortex lattices 
in conventional superconductors, apart from the presence 
of the susceptibility \(q) which reduces to a constant at 
small wavevector and therefore should not qualitatively 
alter the long-wavelength behavior of the vortex lattice. 

Before closing this section, we note that in the expres- 
sion of the elastic Hamiltonian of the SV lattice, (Eq. 
(2.26)), $Q,^(q) for a triangular lattice can be written in 
the form: 



3^/3 (q) 



$ L (q)P Q i /3 (q) + a> T (q)Pj /3 (q), 



(2.50) 



where P^j(q) = q a qp/q 2 ± and Pj^q) = S afj - P^(q) 
are the longitudinal and transverse projection operators 
onto the direction of q^ = q^ /q±_ and the perpendicular 
direction in the (q x , q y ) plane respectively, and where the 
longitudinal and transverse parts of the elastic matrix are 
given by 



*i(q) = cn(q) q]_ + c 44 (q) q 
$r(q) = ceeql +c 44 (q) q 2 . 



2 

z ' 



(2.51a) 
(2.51b) 



Having established the harmonic elastic theory for the 
spontaneous vortex lattice, we next use it to study the 
effects of thermal fluctuations and quenched disorder. As 
discussed in the Introduction, we anticipate these to be 
significantly enhanced relative to that of a conventional 
vortex lattice, due to the softness (vanishing of c 44 ) of 
the spontaneous vortex solid. 



III. THERMAL FLUCTUATIONS AND MELTING 
OF A "SOFT" SPONTANEOUS VORTEX 
LATTICE 

As usual the effect of thermal fluctuations on vortex 
lattice translational order is encoded in the vortex line 
mean-squared fluctuation (u 2 ) relative to the ideal lattice 
positions. This quantity determines the level of transla- 
tional order, namely the intensity of the Bragg peaks (in 
e.g., neutron scattering) through the Debye- Waller fac- 
tor: 



-Q 2 (u 2 



(3.1) 



with Q the length of the shortest reciprocal lattice vector. 

Similarly the location of the melting transition B m (T) 
in the B-T phase diagram can be roughly located by us- 
ing (u 2 ) through the Lindemann criterion 



corresponding to thermal root-mean squared fluctuation 
of vortex lines becoming comparable to the average vor- 
tex spacing a — y^tpo/B, with the Lindemann constant 
cl conventionally taken to be on the order of 0.1. 

Within the harmonic approximation (which, as we will 
find a posteriori, remains valid at low T) (u 2 ) can be 
easily computed using Eq. (2.26) above. Equiparti- 
tion (i.e., simple Gaussian integration) gives the stan- 
dard expression 36 (we use units such that Boltzmann's 
constant ks = 1): 



(u 2 ) 



T 



C66<?i + c 44 (q)<7 



2 ' 



(3.3) 



where we have assumed that en > c 6 6 (as is usually 
the case for materials with a high Ginzburg ratio A/£). 
It is clear from Eq. (3.3) that (as in conventional vor- 
tex lattices 36 ) thermal rms vortex fluctuations are domi- 
nated by the largest wavevectors near the Brillouin zone 
boundary, i.e. near q± ~ A = n/a, q z ~ l/£. At these 
short length scales we expect the spontaneous vortex lat- 
tice to be characterized by a finite tilt modulus, i.e., for 
c 44 (?z) to approach the expression given in Eq. (2.46), 
which should not be too different from that of a conven- 
tional (external magnetic field-induced) vortex lattice. 
We therefore conclude that thermal vortex fluctuations 
(u 2 ) and all other physical observables, such as the melt- 
ing curve and the strength of Bragg peaks, to be quali- 
tatively (and roughly even quantitatively) unchanged by 
the long-wavelength "softness" (c 44 (g 2 — > 0) rts nq 2 ) of 
the spontaneous vortex lattice. We can compare the SV 
lattice fluctuations with the conventional vortex (CV) 
lattice fluctuations ((u 2 )cv ) by treating \(q) pertur- 
batively (valid deep in the ferromagnetic phase where 
|a| » 1) 



(u 2 ) ~ (u 2 ) cv + 



rp CV 

1 c 44 



(q)? 2 



4?rx(q) 



q (c 6 6<zi+c 4 7 4 y (q)^) 2 l + AV 



(3.4) 



(u 2 ) 



.2 „2 



Crfl 



(3.2) 



where is the tilt modulus for a convention vortex 
lattice. This shows that to lowest order in x(q) the SV 
lattice fluctuations are larger, as we would expect due to 
the additional softness. 



IV. RANDOM PINNING WITHIN THE 
HARMONIC ELASTIC THEORY 

As in a conventional vortex lattice, random pinning 
couples to the vortex density p(r) (positional disorder) 
and to the local vortex tilt-field d z u (tilt-disorder, which 
we discuss in more detail in Sec. IV C), that are propor- 
tional to the z- and _L-components, respectively, of the 
local magnetic induction vector field B. The effects of 
disorder on the spontaneous vortex lattice can therefore 
be incorporated through the Hamiltonian H = H e i+Hd, 
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where the elastic energy is given by Eq. (2.26) and the 
pinning energy is given by ' 



I 



dr 



V d (r)p(r) + h(r) • 5 z u(r) 



(4.1) 



where V d (r) and h(r) are the random pinning potentials 
for positional and tilt disorders respectively. 

As we will see in a moment, unlike the effects of ther- 
mal fluctuations studied above, quenched disorder leads 
to vortex lattice distortions that are dominated by long 
length scale fluctuations. As a result, we anticipate that 
the softness associated with the vanishing of the long 
wavelength tilt modulus in a SV lattice will have qual- 
itatively important effects. Hence, to study pinning ef- 
fects, that will be the focus of the rest of the manuscript, 
we specialize to the long wavelength form of the elastic 
Hamiltonian. This is characterized by the elastic ma- 
trix $ Q(3 (q) = $ L (q)Pf (3 (q) + $ T (q)Pj /3 (q) with the fol- 
lowing long wavelength behavior of the longitudinal and 
transverse components: 



*i(q) = c n q± + Kq 4 z , 
$ T (q) = c m q\ + nq* . 



(4.2a) 
(4.2b) 



This elastic matrix is unique to the spontaneous vor- 
tex lattice, with the curvature modulus k given by Eq. 
(2.45b). 



A. Perturbative treatment of random pinning: 
Larkin approximation 

The simplest (yet quite revealing) approximation in a 
study of random pinning of a periodic medium (e.g., a 
vortex lattice) is the Larkin approximation. This ap- 
proximation amounts to Taylor-expanding the positional 
pinning energy (which is a nonlinear function of u) to 
linear order in the phonon displacement u. This can be 
done by expressing the vortex density p(r) as an explicit 
function of u 



p{v) =Rc^p Q e lQ ( r+u ( r », 
Q 

which leads to a positional pinning energy: 
H dp = J drRe^t/ Q (r)e iC > u M , 



(4.3) 



(4.4) 



where the Q's are the reciprocal lattice vectors of the 
vortex lattice and 



U Q (r)=p Q V(r)e^ . 



(4.5) 



We expect correlations of V d (r) to be short-ranged. Even 
if they are not, those of Uq(t) are, due to the e iC *' r fac- 
tor, as discussed in Ref. 29 . This means that one can 
accurately capture the long-length scale physics of the 



problem by taking them to be zero-ranged and Gaussian, 
characterized by a second cumulant: 



[/ Q (r)^,(r')=A Q ^(r-r')5 Q , Q - 



(4.6) 



where the set of phenomenological parameters Aq's de- 
pend upon the microscopic natures of the impurities, the 
vortex lattice and their interactions. 

In using the Larkin approximation we expand the peri- 
odic nonlinear piece for small u, valid only at short length 
scales such that the typical induced vortex displacement 
u ms = (u 2 ) is small compared to the correlation length 
of the random potential. With this approximation we 
can write the pinning energy as: 

H d w J dv [f (r) • u(r) + h(r) • z ii(r)] , (4.7) 

with f(r) an effective random force with two-point cor- 
relations: 



f a {v)f ti (v')=A f 5 ati 8{v-v'). 



(4.8) 



where A/ = J2q Q 2 Aq. This leads to a total Hamilto- 
nian in this short-scale Larkin regime that is given by: 



H = H el + Jdr f(r)-u(r). 



(4.9) 



Above we have dropped the tilt disorder h(r), since it 
is clearly subdominant to the uncorrelated random force 
f(r) arising from the Larkin approximation to the posi- 
tional pinning disorder. Standard Gaussian integration 
then leads to the disorder-averaged phonon correlation 
function 

C"(r) = C{{t) + C£ (r) = <[u(r) - u(O)P), (4.10) 

where the longitudinal and transverse parts (r) (with 
s = L, T labeling the longitudinal and transverse polar- 
izations, respectively) of (r) are given by: 



C/(r) = / 



2A, 



2T[1 - cos(q ■ r)] 
1 — cos(q • r) 



L 



[B s ql 



nq z 



l2 ' 



2A/[1 -cos(q- r)] 
q [B s ql + m i] 2 

(4.11) 



with Bl = Cn and Bt = cqq the bulk and shear moduli 
of the vortex lattice, respectively, and where, in going 
from the first to the second equation, we have neglected 
the thermal contribution to C/(r), that is subdominant 
to the disorder part at small wavevectors. 

Simple dimensional analysis shows that, for d < 9/2, 
the above integrals are dominated by long length scales, 
justifying our use of long-wavelength elastic theory de- 
rived in the previous section. In d dimensions we find: 



C/(r = R ± ) 
Cf(r = R z z) 



A/ 



K (5-d)/2 B (d-l)/2 



9-2d 



(4.12a) 
(4.12b) 
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where we have dropped overall numerical factors of order 
unity, and z and _L denote directions along and perpen- 
dicular to the flux lines, respectively. In the physically 
relevant case of d = 3, the above two expressions be- 
come: 



C{{R±) 
Cf{R z ) 



A f R 



3/2 



A f Rl 



(4.13a) 
(4.13b) 



Equating these correlation functions to the square of the 
lattice spacing, a 2 , leads to highly anisotropic Larkin 
lengths in the (xy) and z directions, given respectively 
by: 



R„ 



R z . 



V A, ) ' 



A/ 

' K B,a 2 \ i/3 



(4.14a) 
(4.14b) 



These lengths characterize the dimensions of ordered 
Larkin domains beyond which pinning dominates over 
elastic energy, disrupting the translational order of the 
spontaneous vortex lattice. The finiteness of these 3d 
Larkin domains demonstrates the absence of long-range 
translational order even for arbitrarily weak disorder. By 
this last criterion, the Larkin lengths also define the range 
of length scales over which the above random-force per- 
turbation theory, as defined by the Hamiltonian (4.9), is 
valid. At longer length scales the rms value of u exceeds 
the lattice spacing, which invalidates the approximation 
we used in obtaining Eq. (4.7). Also, because the ran- 
dom force dominates over the random tilt disorder, it 
is the positional pinning (rather than tilt disorder) that 
determines the size of Larkin domains. 

However, as is well-known for conventional vortex lat- 
tices and other periodic media 38 , on longer length scales 
the Larkin approximation highly overestimates the effect 
of disorder. On scales beyond Rjr and R z c , a more sophis- 
ticated approach that takes into account the nonlinearity 
of the positional disorder is necessary. Next we shall use 
the replica Gaussian variational method 39,34 to treat the 
long-scale effects of the (seemingly dominant) positional 
disorder alone. For now, ignoring tilt disorder and non- 
linear elasticity, we shall find, that, in spite of the soft 
elasticity of the tilt modes in a spontaneous vortex lat- 
tice, positional disorder leads to a logarithmic growth of 
the displacement correlation function (r) in the phys- 
ical case of d = 3 dimensions, as it does in conventional 
vortex lattices in ordinary superconductors. 38,34 



B. Positional pinning on long length scales 

As discussed above, the perturbative (random-force) 
treatment of positional disorder breaks down on length 



scales longer than the Larkin lengths R^r and R z c of Eqs. 
(4. 14a)- (4. 14b). In this Subsection we shall use a replica 
variational analysis 40,41 similar to the one carried out in 
Refs. 39, 34 to find the contribution to the long distance 
behavior of the two-point correlation function for rela- 
tive vortex displacement from the seemingly dominant 
positional disorder alone. Our starting point here is the 
Hamiltonian: 



# = jT^[tfgi + /«£]|u(q)| 

+ J drReJ2UQ(r)e lQut 



(4.15) 



where for simplicity we have used isotropic elasticity with 
a single in-plane elastic constant K and the positional 
disorder correlations are given by Eq. (4.6). Employ- 
ing the standard replica "trick" 42 to average over the 
quenched disorder 38 , we obtain the following replicated 
effective Hamiltonian 

H *ff - E / \ [ K & + K «* ] u <^) • u ^- f i) 

m f A 

"EE/ dr ^F cos[Q.(u (r)-u fe (r))], 



a, b=l Q#0 l 



(4.16) 



with a, b labeling the m — > replicas, Q denoting recip- 
rocal lattice vectors, and Aq defined by Eq. (4.6). In the 
above we have also dropped the coupling of disorder to 
the long-scale fluctuations in the vortex density — /o V-u, 
that in 3d is subdominant to the short-scale pinning that 
we have kept. 

We now want to study the above nonlinear Hamilto- 
nian via a Gaussian variational approximation 39,34 with 
the following harmonic trial Hamiltonian: 



H„ — 



/(G _1 (q)Wua(q)-u 6 (-q), (4.17) 

Jq 



which is parameterized by a variationally-determined in- 
verse "propagator" 



[G- 1 (q)] ab = (Kql + K q 4 z )5 ab -<J ab . 



(4.18) 



characterized by an m x m self-energy matrix a a b of vari- 
ational parameters, that encodes the average vortex po- 
sitional correlations. 

Minimization of the variational free energy 



(H, 



eff 



Hi, 



T In Z„ 



(4.19) 



where Z v = Tr( cxp(— H v /T)) and (■■■)„ denotes aver- 
aging with statistical weight exp(—H v /T)/Z v , leads to 
saddle-point equations that are very similar to those de- 
rived in Ref. 34. As discussed in more detail there, the 
long-scale properties of the pinned state are characterized 
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by a replica-symmetry broken matrix A a b with "hierar- 
chical" structure, that in the m — > limit can be viewed 
as a single function o~(v) of a real variable v in the inter- 
val [0, 1], that is the key quantity of the replica method. 
The hierarchical matrix a(v) is a solution of the saddle- 
point equation (here Qo is the magnitude of the smallest 
reciprocal lattice vectors): 



a(v) f 

Jq 



TQl 



[Kqi + Kql + [a]{v)\ 



= 1. 



(4.20) 



where we denote by [c](v) the quantity [cr](w) = vo~(v) — 
J " du <t(u). The solution of the saddle-point Eq. (4.20) 
proceeds in much the same way as for a conventional, 
field-induced vortex lattice. 34 We refer the reader to Ap- 
pendix A where we give more details on this solution, 
and present here only the final results for the disorder- 
averaged correlation function C p (r) = ([u(r) — u(0)] 2 ) 
(we use the superscript p in C p (r) to distinguish the long- 
scale translational correlator in presence of the periodic 
density pinning alone - i.e., in the absence of tilt disorder 
and clastic nonlincaritics from the quantity (r) of the 
previous paragraph which corresponds to the Larkin ap- 
proximation of random pinning forces and which is only 
valid inside the Larkin domains defined by the length 
scales and R z c of Eqs. (4.14a)-(4.14b)). 

For a spontaneous vortex lattice in three spatial di- 
mensions, we find (see Appendix A): 



C. Random tilt disorder 



In addition to the usual effects of pinning leading to a 
random coupling to the vortex density, Eq. (4.15) (corre- 
sponding to the pinning of the ^-component of the mag- 
netic field), it is essential, in studying a spontaneous vor- 
tex solid, to include a random tilt disorder, as appearing 
in H4, Eq. (4.7). As schematically illustrated in Fig. 
2, physically, such a random term corresponds to local 
random torques being exerted on the vortex lattice by 
random, short-ranged correlated clusters of pinning sites. 
The interaction with these randomly oriented clusters of 
point pinners is similar to the interaction of vortex lines 
with a randomly oriented local transverse magnetic field 
h(r) (not to be confused with the screening field around 
vortices of Sec. II), that couples to the transverse com- 
ponent of the magnetic induction vector field Bj_, and 
therefore leads to random coupling: 



Hi 



-! 



dr h(r) • <9 z u(r) . 



(4.23) 



For simplicity we take h(r) to be a Gaussian random field 
with zero mean and correlations: 



Mr)M r ') = A t <J Q/3 <5(r - r') . 



(4.24) 



ln(Ar_L) . 



(4.21a) 
(4.21b) 



where (here d± = d — 1 and Cd ± is a numerical constant 
defined in Appendix A) 



u 



3V2c d± TQ 2 

4(7-2^)^/4^3/2 • 



(4.22) 



We thus see that the softness of the tilt modes of the 
spontaneous vortex lattice does not affect the large dis- 
tance behavior of the correlation function C p (r), which, 
like in ordinary (field-induced) disordered vortex lattices, 
grows only logarithmically at long length scales. This 
is expected on general grounds, since (weak) disorder- 
induced logarithmic growth is a consequence of the pe- 
riodic form of positional pinning in Eq. (4.16), which 
is guaranteed by the identity-symmetry of vortex lines. 
However, this result is valid only if we ignore the tilt 
disorder and nonlinear elasticity, which, as we will show 
below become important in d = 3 < 7/2 dimensions. 
In what follows, we shall consider the effect of random 
tilt disorder on the spontaneous vortex lattice, and show 
that it leads to power-law growth of the vortex line dis- 
tortions in d = 3 dimensions, and therefore, at long scales 
strongly dominates over the logarithmic growth induced 
by positional disorder found above. 



Independent of the microscopic mechanism, it is simple 
to see, that tilt disorder is always generated from the po- 
sitional disorder upon coarse-graining even if it is left out 
in the original model. Such disorder is generally left out 
in a treatment of conventional vortex lattice, as in three 
dimensions it (as well as the aforementioned coupling of 
disorder to the long-scale fluctuations in the vortex den- 
sity — V(r)p V • u) is subdominant to the positional pin- 
ning that we have studied in the previous Subsection. 




FIG. 2. Schematic representation of the mechanism of ran- 
dom tilt pinning. Randomly oriented anisotropic clusters of 
point-pins exert random torques on vortex lines. 
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However, as we will show below, because of the softness 
of the spontaneous vortex lattice (more specifically, the 
vanishing of its long wavelength tilt modulus C44), this 
tilt disorder leads to power-law distortions of the sponta- 
neous vortex lattice in three dimensions. It is therefore 
the dominant form of disorder, that must be taken into 
account on long scales. 

To see this we treat tilt disorder alone. Because tilt 
disorder couples only linearly to the phonon u field (in 
the absence of positional disorder and nonlinear elastic- 
ity) it is straightforward to study its effects on the vor- 
tex positional correlations. Using the Hamiltonian H = 
H e i+H t , a standard calculation gives for the longitudinal 
and transverse co ntributions of th e positional correlation 
function C*(r) = ([u(r) - u(0)] 2 ) = C£(r)+C^(r). Here 
we use the superscript t to distinguish this correlator as- 
sociated with tilt disorder from the correlation functions 
and C p of the previous two Subsections that were as- 
sociated with positional disorder. Keeping in mind that 
the subscript s denotes longitudinal (L) and transverse 
(T) contributions we find that: 



Cl(r) = 2T 



+ 2A 



2A ; 



I 

L 



q\ [l - cos(q ■ r)] 
[nq* + B s q\] 
q 2 z [l - cos(q ■ r)] 
q [nqi + B s qlY 
q 2 z [l - cos(q • r)] 



q [ K qi + B s qlY 



(4.25) 



where, in going from the first to the second equation, we 
have dropped the subdominant thermal part and only 
kept the T — random pinning contribution. Asymp- 
totic analysis of the above integrals gives (we again drop 
overall numerical factors of order unity): 



C2(r±) 



Cl(z) 



A t rj /2 

«3/4iJ«/4 ' 

A t 



.,7-2d 



K (5-d)/2 B (d-l)/2 
A t Z 



kB s 



(4.26a) 

(4.26b) 
(4.26c) 
(4.26d) 



where in Eqs. (4.26b) and (4.26d) we have specialized to 
the physically relevant case of three dimensions. Defin- 
ing tilt-disorder spatial coherence lengths as the 
lengths for which C*^,^ 2 ) = a 2 , we find 



B 5 J\^ 



f nB s 



(4.27a) 
(4.27b) 



We stress that (in the absence of positional disorder, elas- 
tic nonlincarities, and dislocations) the power-law growth 



of the correlation function Cj(r) (Eqs. (4.26b)-(4.26d)) 
extends out to arbitrarily long scales. This is in contrast 
to the Larkin random force approximation of Sec. IV A. 
Hence at long length scales, for d = 3 < 7/2, C*(r) quite 
clearly dominates over the much slower, logarithmically 
growing distortions C p (r) (Eqs. (4.21a)-(4.21b)) created 
by translational disorder alone. This is a very important 
point and is one that cannot be seen just by comparing 
the coherence lengths (£t) in the presence of tilt disor- 
der alone (given by Eqs. (4.27a) and (4.27b)) with those 
(£ p ) in the presence of positional disorder alone (given 
by Eqs. (4.14a) and (4.14b)). For equally weak disor- 
der strengths (i.e. small A t and A/) one can see that 
£i 3> £ p and Larkin length is determined by the posi- 
tional disorder, i.e. £ p = R c . However, at length scales 
beyond this Larkin length (where the Larkin approxi- 
mation breaks down) the vortex lattice distortions are 
controlled by the random tilt (rather than positional) 
disorder. 43 Therefore, henceforth we shall ignore transla- 
tional disorder altogether and focus on tilt disorder only. 
As we will show below, despite this considerable simpli- 
fication, the asymptotic properties of the spontaneous 
vortex lattice remain nontrivial and rich due to the in- 
terplay of tilt disorder and anharmonic elasticity that we 
have ignored so far. 



V. PERTURB ATI VE TREATMENT OF ELASTIC 
NONLINEARITIES AND TILT DISORDER 

The elasticity of a pinned line-crystal can be most eas- 
ily formulated in terms of the (Lagrangian) left Cauchy- 



Green strain tensor v, 



'a/3, 



defined by 



Va/3 



-(diR a diRp - 8 a p). 



(5.1) 



For a line-vortex solid characterized by a conformation 
R(X M ,z), defined in Eqs. (2.8), (2.24), the nonlinear 
Lagrangian strain tensor reduces to 



V a f3 



5(* 



up + dfju a + d l u a d i up 



(5.2) 



with a, (3 ranging over the (i^-dimensional subspace (a;, y) 
transverse to the average spontaneous vortex line direc- 
tion, which we take to be z, and i ranging over the full 
space (x, y, z). In conventional elastic solids, for weak dis- 
order and low temperature the gradients of the phonon 
field u(r) are usually small and it suffices to approxi- 
mate the full nonlinear strain tensor by its harmonic 
part = \{d a up + dpu a ). However, there exists a 
novel class of "soft" elastic systems 22 in which nonlinear 
elasticity plays an essential qualitative role. A unifying 
feature of solids in this class is their underlying, spon- 
taneously broken rotational invariance, that strictly en- 
forces a particular "softness" of the corresponding Gold- 
stone mode Hamiltonian. As a consequence, the usually 
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small nonlinear clastic terms are in fact comparable to 
harmonic ones, and therefore must be taken into account. 
Similar to their effects near continuous phase transitions, 
but extending throughout an ordered phase, fluctuations 
and disorder drive elastic nonlinearities to qualitatively 
modify such soft states. The resulting strongly inter- 
acting ordered states at long length scales exhibit rich 
phenomenology such as a universally nonlocal elasticity, 
a strictly nonlinear response to an arbitrarily weak per- 
turbation and a universal ratio of wavevector-dependcnt 
singular elastic moduli, all controlled by nontrivial in- 
frared stable fixed points, illustrated for the problem at 
hand in Fig. 3. As we have recently argued in a brief 
communication 30 a spontaneous vortex solid, with its 
symmetry-enforced vanishing tilt modulus C44 is a mem- 
ber of the class of such soft solids. Hence, as we will show 
below, to understand its pinned state the nonlinear part 
of the strain tensor v a p, Eq. (5.2) must be taken into ac- 
count. In this Section we demonstrate the importance of 
nonlinear elastic terms via a simple perturbation theory 
in these nonlinearities and tilt disorder. 

The full nonlinear elastic Hamiltonian of a spontaneous 
vortex lattice pinned by the tilt disorder is given by (here 
the Lame elastic constant A should not be confused with 
the London penetration depth): 



is the standard harmonic elastic part, while 



H 



dr 



1 



A 



-k (d z z u) z + - v z aa + h(r) • 9 z u(r) 



(5.3) 



where k is the curvature modulus, and fi = c 6 6 and 
A = Cn — c 6 6 are Lame coefficients, respectively char- 
acterizing the in-plane shear and in-plane bulk moduli 
of the spontaneous vortex solid. The key feature that 
distinguishes the spontaneous vortex solid from a conven- 
tional field-induced one, and requires us to keep nonlinear 
terms in d z u a , is the absence of the n z ±v za out-of-planc 
shear term. This is because of the exact vanishing (en- 
forced by rotational-symmetry) of the tilt-shear modulus 
/j, z ± = C44. In contrast, because /i and A arc finite, the 
corresponding elastic nonlinearities in d a up need not be 
retained, as they are clearly subdominant at long length 
scales. Consequently, for a spontaneous vortex solid the 
nonlinear elastic strain that appears in iJ, Eq. (5.3), is 
given by: 

Va/3 ^ ^ (d a Ufj + d (i U a + dzUadzltpj . (5.4) 

The Hamiltonian H above can be written as: 



H = H + Hi + Hd , 



(5.5) 



where 



«,-/ 



dr 



\n{d 2 z n) 2 + \n{d a upf 



+ \ ii{d a u fj ){dpu a ) + i \(d a u 



lid a u (i d 7 u a d z u (} + ^ {d a u a ){d z Uf)f 



+ -(A + 2/i) {d z u a f{d z upf 
is the nonlinear elastic contribution, and 
H d = Jdr h(r) • d z u(r) 



(5.7) 



(5.8) 



is the previously-introduced random tilt disorder. The 
quenched tilt disorder can be most efficiently treated via 
the standard replica trick, 42,46 



Z n 

Z n = lim 



n^O n 



(5.9) 



which allows us to formally average over disorder and 
therefore work with a more convenient translationally in- 
variant field theory. Introducing n replica fields labeled 
by a and averaging the resulting "partition function" , Z n 
over the tilt disorder, h(r), we find 



Z n = J [dui(r)] • • • [du n {r)] C -P( H ^+ H ^) j 
where /3=1/T, and 

n 



(5.10) 



(5.11) 



o=l 



and we denote by Ho n the following (quadratic) Hamil- 
tonian: 

H 0n = Yl [ ^(q)<(q)^(-q), (5.12) 

a,b q 

with the kernel: 

-\a& 



A 



r$(q) = $ a/S (q)J„6 - =f q 2 z 5 a p 



(5.13) 



where $ Q/3 (q) = $ L (q)Pf /3 (q) + $ T (q)Pj /3 (q) and $ L (q) 
and < I ) T(q) are given by Eqs. (4.2a) and (4.2b). From 
the above equation the noninteracting elastic propaga- 
tor G a p(q) = tt^ 1 (u a (q)uf3(— q))o (O being the system 
volume) can be easily obtained (see Appendix II) and is 
given in the limit n — > by: 



0<l) ^ Gf (q)P^(q) + G5?(q)P^(q) , (5.14) 
where 

Gf (q) - 2T7/(q) 5 ab + A t q 2 z Tf( q ) , (5.15a) 
G a T b (<D = TT^q) S ab + A t q 2 z r T 2 (q) , (5.15b) 

(5.6) 

with the kernels: 
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TL(q) = (A + 2[i)q\ + nq 

iMq) = mil + m\ ■ 



Z 1 



(5.16a) 
(5.16b) 



To assess the importance of elastic nonlinearities we 
coarse-grain Z n . To this end we decompose the displace- 
ment field u(r) into low and high wavevector parts: 



u(r) = u<(r) + u>(r) . 



(5.17) 



We then perturbatively (in the nonlinear elastic terms) 
integrate out the quickly varying, short-scale components 
u > (r), and thereby obtain a representation of Z n : 



Z* = J [du< (r)]c-^„ Z > (c-P H ^ 
= j[du<{v)\ e~^, 
in terms of a coarse-grained Hamiltonian H e g 
H cfi [u] = H<Ju] + (H ltn ) Q> - ^{Hl n ) c 0> 



o> 



(5.18) 



(5.19) 



where Zq = Tr(— /Jifgn), (• • -)o> indicates an average 
with statistical weight exp(—/3HQ n )/ZQ, (i?i„)o> de- 
notes the connected cumulant = (H^ n ) 0> — (if^) 2 ^, 
and where short-scale modes extend out to wavelengths 
of size L z x L±. 

Restricting our attention here to the leading order 29 
nonlinear part of H ln : 



H,, 



a=l 



dr {iid a u%d z u a J) z u% + - d a u a a {d z u%) 



(5.20) 

a standard calculation (with details relegated to Ap- 
pendix III) shows 47 that the resulting effective Hamil- 
tonian, i/ c ff[u < ], can be put into the same form as the 
original Hamiltonian, i?on[u], but with effective elastic 
moduli perturbatively corrected by the coarse-graining 
procedure. For example, as we show in great detail in 
Appendix III, the lowest order correction (5/x to the bare 
clastic shear modulus \i is given by: 



5jJL ■ 



Td ± (d± + 2) 



ft J 



d d ~ 



q.L 4 

Qz 



(2ir) d 



+ (d\ - 2)(G a T b (<DY + 2d ± Gf (q)GS?(q) 



2(Gf(q)) 2 
(5.21) 



where G°£ and G^ are the noninteracting propagators of 
Eqs. (5.15a)-(5.15b) and where the ">" sign on the q z 
integral indicates that we impose the IR cutoff 1/L Z1 re- 
stricting the integration range to the region \q z \ > 1/L Z . 
Evaluating the above expression to leading order in the 
IR cutoff L z (and ignoring the subdominant thermal 
parts) gives the perturbative relative corrections to Lame 
elastic moduli: 



ji A 



'7 -2d 



(5.22) 



that for d < 7/2 diverge with L Z ,L±. As we will see in 
subsequent Sections (with details given in Appendix III) 
similar analysis shows that the perturbative corrections 
to the curvature modulus k and the tilt-disorder strength 
A t also diverge with system size like L 7 z ~ 2d . Thus, at long 
length scales perturbation theory breaks down for spatial 
dimensions smaller than the upper critical dimension 



due — ~ • 



(5.23) 



Note that this is the same as the critical dimension be- 
low which positional fluctuations diverge. This is not a 
coincidence. Rather, it is a consequence of the fact that, 
in power-counting terms, the leading cubic and quartic 
anharmonicities that we have kept have the same num- 
ber of derivatives as the leading harmonic terms. In doing 
this power counting, note that each _L derivative counts 
as two z-derivatives, as can be seen from the strongly 
anisotropic form of the harmonic terms. Given this, the 
only difference in scaling between the anharmonic and 
the harmonic terms must come from the extra powers of 
u in the anharmonic terms. Hence, the anharmonic terms 
only become important when u scales up with increas- 
ing distance; i.e., when u fluctuations grow with length 
scale, as they do below d = d uc = | . 

We define a nonlinear (crossover) length scale at which 
perturbation theory breaks down by the value of L z at 
which the relative corrections to the elastic moduli and 
disorder variance, SX/X, Sfi/fi, 5k/k and 8A t /A t , be- 
come of order 1. Equation (5.22) leads to the following 
estimate for this nonlinear crossover length : 



7-d d-3 



(B 



for d = 3. 



(5.24a) 
(5.24b) 



A corresponding perpendicular nonlinear crossover 
length £± L can also be defined by imposing a 1/Lj_ in- 
frared cut-off on the q_i_ integrals in perturbation theory, 
and is given by 



C NL 



K / 



(5.25) 



Hence, we find from the above perturbative coarse- 
graining procedure that at sufficiently long length scales 
(greater than mm[(^ L ,^ L )]), the anharmonic elastic- 
ity becomes qualitatively important, even for arbitrar- 
ily weak pinning, thus invalidating the description of 
the pinned spontaneous vortex solid by a conventional 
harmonic elasticity theory. In order to describe the 
pinned spontaneous vortex solid on scales longer than 
(£± L ,£, z srL ), we need to elevate the above perturbation 
theory to a renormalization group analysis. We turn to 
this task in the next Section. 
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VI. NONLINEAR ELASTICITY: RG ANALYSIS 

Having established the importance of the nonlin- 
ear elasticity (as indicated by the divergent perturba- 
tion theory of previous Section), we now employ the 
standard Wilson momentum-shell renormalization-group 
method 48 to study the pinned spontaneous vortex solid 
on length scales longer than £nl- The results we derive 
here were first reported in Rcf. 30. 



A. Recursion relations and zero-temperature fixed 
points 

To this end, we use phonon mode decomposition, Eq. 
(5.17), with high wavevector modes u> restricted to an 
infinitesimal cylindrical shell: 

A z e- £ < \q z \<A z , (6.1a) 
0<|gj<oo, (6.1b) 

where £ is taken to be an infinitesimal positive number. 
The procedure is to integrate these short length scale 
modes u > (r) out of the replicated partition function Z n 
and to interpret the result in terms of a Hamiltonian of 
the same form but with effective ^-dependent parame- 
ters. To simplify the algebra that follows we introduce 
anisotropic rescalings of lengths and fields: 

(6.2a) 
(6.2b) 
(6.2c) 



z = e e z' , 
u<(r)=e^<(r'). 



This leads to zero-th order (in elastic nonlincaritics) re- 
cursion relations 



M (^ M c« d - 3 ^+ 1+2 ^, 
— (£) = ^ e ((d ~ 1)w ~ 1+20)£ 



(6.3a) 
(6.3b) 
(6.3c) 

j j (6.3d) 

It is convenient (but not necessary) to take advantage 
of the rcscaling freedom (6.2a)-(6.2c), and to choose the 
arbitrary exponent </> to be 

<j> = 2-u>, (6.4) 

so as to keep fixed 49 the form of the nonlinear strain 
tensor v a/3 = (d a u f3 + dpu a + d z u a d z up)/2. 

Integrating out the high wavevector modes u> to sec- 
ond order 47 in elastic nonlinearities leads to the follow- 
ing corrections to effective Hamiltonian parameters (with 
calculation details provided in Appendix IIIB): 

(6.5a) 
(6.5b) 
(6.5c) 

(6.5d) 



SD = gF D (x)D de, 
Sk = gF K (x)n dl , 



de. 



where D = A + fi and we have defined dimcnsionless cou- 
plings x and g such that: 



A 

x = — . 
i 1 



9 = 



\/ 2Crf-i 
"64 



d — 7 3-d 



A t k 2 fx 2 A ; 



2d-7 



(6.6a) 
(6.6b) 



the latter controlling the perturbative expansion; and 
where the scaling functions Fi(x) are given by: 



F D (x) 



+ 



30(x+ l) 2 
90s 3 + 234a; 2 + 144a; - 64 



135a; 3 + 351a: 2 + 316a; + 164 + 



(2 + a;) 5 /4 



F»{x) = 
F K (x) = 



15(x+ 1) 
1 



97 + 17 x - 



5(x + l) 
1 

10(x + 1) 



1686 + 366a;- 



(6.7a) 

(x + 2f/^ ^ 
3896+3156a;+580a; 2 i 



152 + 72a; 



(x + 2) 5 / 4 



(6.7c) 



543 + 159a; 



614 + 230a; 



(x + 2) 1 / 4 



. (6.7d) 



Combining these with the zero-th order rescalings, Eqs. 
(6.3a)-(6.3d), and Eq. (6.4), leads to the differential flow 
equations: 



dD 

H 

d\i 

~d~e 

dn 

dl 
dA t 

de 



It is easy to verify that the above differential flow equa- 
tions, Eqs. (6.8a)-(6.8d), lead to the following two closed 
flow equations for the two dimcnsionless couplings g(i) 
and x(e): 



((d- 


5)u + 5 + gF x (x))D(£), 


(6.8a) 


((d- 


5V + 5 + 5 F M (x)) M (£), 


(6.8b) 


((d- 


3)uj + 1 + gF K (x)) n(e) , 


(6.8c) 


((d- 


3)uj + 3 + gF At (x))A t (e). 


(6.8d) 



dg_ 

de 



(7-2d)+g(F At (x) + ^F,(x) 



+ d -l FK (x)) 



9- 



(6.9a) 



^=g[F D (x)-F,(x)](x(e) + l). (6.9b) 



The fixed point of Eq. (6.9b) (i.e., such that (dx/d£) = 
0) can be easily computed by noting that for finite g they 
are given by the zeros of the function (F\(x) — F^(x)) 
plotted in Fig. 4. These are given by: 
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-4/5 



-0.033 



FIG. 3. Renormalization group flow of the ratio of the 
Lame coefficients x = X/fi. The fixed point at x 2 = —0.033 
controls the long-scale properties (anomalous elasticity and 
associated phenomenology) of the pinned spontaneous vortex 
solid. The critical point at x 1 = —4/5 controls the mechanical 
instability of the spontaneous vortex solid. 



Xl 
X2 



-4/5, 
-0.03272. 



(6.10a) 
(6.10b) 



As can be seen from Fig. 4, (Fd(x) — F^(x)) is positive 
for xi < x < X2, and negative otherwise. Therefore x\ = 
—4/5 is an unstable fixed point, while x^ — —0.03272 is 
the stable fixed point of interest to us (see Fig. 3). 

At the stable X2 — —0.03272 fixed point, that we shall 
denote hereafter by a:*, the recursion relation (6.9a) for 
the dimensionless coupling g(i) reduces to 



^= [2e- 12.8359 5 ] ff , 



(6.11) 



with e = 7/2 — d controlling the validity of the expansion. 
This implies a stable fixed point (to 0(e 2 )) 



(a;*, 5*) ~ (-0.03272, 0.1558 e) . 



(6.12) 



that controls long-wavelength, weak-disorder and low- 
temperature properties of the pinned spontaneous vortex 
solid. 

The existence of this nontrivial fixed point has dra- 
matic consequences for the long-scale elastic properties 
of the pinned SV lattice. As discussed above, it leads to 
the anomalous elasticity (common to soft solids), that we 
explore below. 

The critical point x\ = —4/5 and its instability cor- 
respond to the stability limit, in 7/2 dimensions, of the 
elastic constants, A and /i, beyond which the vortex lat- 
tice is unstable against bulk distortions, i.e. u(r) cx Te- 
stability of the vortex lattice against bulk distortions re- 
quires that the bulk modulus, B = A + (2fi/d*), (not to 
be confused with magnetic flux; d* = d — 1 is the di- 
mensionality of the vortex lattice) be positive. At the 
stability limit, i.e. when B = 0, 



A 

x = — 



d- 1 



(6.13) 




FIG. 4. A plot of the function (F D (x) - F M (x)), with its 

zeros, (Fd(x) - F M (x)) = at xi = -4/5 and x 2 ~ -0.03272 
determining the fixed point solutions for the Lame coefficients 
ratio, x* = A* / fj,* . 

which for d = 7/2 gives x = —4/5, the value of the 
unstable fixed point. That the recursion relations re- 
flect this stability limit is a nontrivial check that the A 
and fj, graphical corrections were correctly evaluated. For 
x < —4/5 the flows run away negatively, corresponding 
to a phase transition from the stable vortex lattice phase 
to one which, presumably, has different lattice constants. 
The associated phenomenology would be an interesting 
subject for future research. 



B. Long length-scale correlation functions: matching 

Many of the basic properties (length scale-dependent 
elastic moduli and the associated anomalous elasticity) 
of the spontaneous vortex glass are encoded in the long- 
scale limit of the 2-point phonon correlation function 

^c^W = (^a b p( c l)) i defined and derived in Appendix 
II for the harmonic elastic case (valid on length scales 
shorter than (,nl)- 

We can utilize the renormalization group transforma- 
tion and the harmonic propagator to compute the cor- 
relation function at longer length scales despite it being 
nonperturbative in the elastic nonlinearities. To this end, 



recalling the definition of G% (q) , 



Gab 
a/3 



(q) 



K(q)^(q')) 
(27r) d 5(q + q') ' 



(6.14) 



together with the rescalings, Eqs. (6.2a)-(6.3d), it is easy 
to see that 



G^(q)=c 



((d-l)w+l+20)£_ 



(<(q ± e^, 92 e^(q' ± e^, 9 y) ) 



(6.15) 



(2Tr) d S(q ± e^ e + q'_ L e^)5(g z e f + q' z e l ) ' 

= e««)^G^(q i e w< ) ^), 

where in going from the first to the second line we made the convenient choice </> = 2 — lo. Finally, integrating out the 
short-scale phonon modes u > , gives a scaling relation satisfied by the correlation function: 

G^(q ± ,q z ,X,^K,g) = S^+^Gf^^ 1 ,q z e l ,\{t),v{t),K{i),g{Z)) . (6.16) 



16 



This scaling (matching) relation demonstrates the real 
power of the renormalization group transformation. 50 
Choosing e l large, we can relate the long-wavelength 
(small q) correlation function (that is impossible to com- 
pute perturbatively because of the appearance of infrared 
divergences in the direct perturbation theory) to the cor- 
relation function at large q, (e.g. at the Brillouin zone 
boundary), that can be easily computed using simple per- 
turbation theory. 



1. Long wavelength behavior of k and A t 

To utilize the above scaling relation we first focus on 
the q± — case and take I* such that 



= A, 



(6.17) 



We then consider sufficiently long length scales, A z /q z 3> 
1, so that l* is large enough for the coupling constants 
((/(£*), x(£*)) to approach their fixed point values (g*,x*). 
Then using Eq. (6.16) together with the zero-th order 
approximation for its right hand side, we find 



G a a b (q± = 0, q z , A, fi, k, g) = e^ d ~ 3 ^ r G^(0, A z , \(t ), n{t ), n(t ), Qt ) 



a ((d-3)w+5)r 



T 



Sab 



Using Eq. (6.17) to eliminate I* in favor of q z , the first term in square brackets gives: 

T rp e ((d-3)oj+5)e* 



0~a(3 ■ 



a ((d-3)w+5)r . 



k{£*)A z 



SabSaf: 



K^(A ;; /9 z ) 4 e«' i - 3 ) u + 1 +9- ,? «( I '))« 
T(A z /q z f d -^ 



SabS, 



a{3, 



mi{A z /q z ) 

T 



4+((d-3)w+l+g.F K (x,)) 



0~abO~af3, 



Kq\(A z jq z y 



(6.18) 

(6.19a) 

(6.19b) 

(6.19c) 
(6.19d) 



We note that as required, the arbitrary rescaling expo- 
nent w dropped out from the physical correlation func- 
tion. In the above expression, we can identify n(q z ) as 
the effective wavevector-dependent curvature modulus: 

«(&) = «(£)~"". ( 6 ' 2 °) 

with the exponent 

Vk= 9*F K {x r ) , (6.21) 

clearly a universal number determined by the fixed point. 
A similar reasoning for the second term in (6.18) leads 
to the wavevector-dependent disorder variance: 

A t (g z )=A t (f;)^ (6-22) 

that corresponds to effective long-range correlated tilt 
disorder, with a universal exponent given by 

VA t =g*F At (x.). (6.23) 

From equations (6.21) and (6.23) and the definitions of 
the functions F K (x) and F& t (x), Eqs. (6.7c)-(6.7d), we 
obtain the following numerical values of the exponents 
r] K and r]A t ■ 



r] K = 1.478e (6.24a) 

~ 0.739 in d — 3 dimensions, (6.24b) 

7] At = 0.414e (6.24c) 

~ 0.207 in d = 3 dimensions. (6.24d) 



2. Long wavelength behavior of A and \i 

Similarly, in order to determine the long wavelength 
behavior of the effective Lame moduli A and fj,, we again 
utilize the scaling relation, Eq. (6.16). Here, quite 
clearly, we need to keep q± finite. However, for con- 
venience, we choose q±£,± L <C {qz£,^ rL )'', so that (as in 
the previous subsection) it is q z that controls the length 
scale dependence of the elastic moduli, and therefore the 
rescaling parameter e e ' is still determined by Eq. (6.17): 

e<- = (!). (6.25) 

The longitudinal part of (q) then satisfies: 

Gt b ( q± ,q z ) = c ^-^+W* Gt(q ± e^,A z ) 
Te ((d-3)oj+5)e* 

= [A(£*) + 2[i(e*)] q 2 j_e 2t * + k(^*)A4 ' 
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(6.26) e ((d-5V+5)r [ Ae s,F,(x,)r + 2 M e^ F ^ a: *) r ] q\e u * 

Utilizing the elastic constants flow equations, Eq. (6.8), Wz ' 

the denominator of the expression on the right-hand side 

of this last equation can be rewritten in the form Using Eq. (6.25) we obtain: 

b T 
G a L (qj_,<fc) = ^ AeS , F , (x , K , +2/ie s.F,(*.)/'] ql + Kq i e g,F K { x ,)i* ' ( 6 - 27 ) 



which shows that at long length scales the longitudinal 
part of the correlation function G£ h (qj_, q z ) takes on its 
harmonic-elasticity form: 



G a L (q_L,g z ;A,^,K,.g) = 



T 



[X(q z ) + 2fi(q z )]ql + K(q z )q^ 

(6.28) 



but with crucial wavevector-dependent elastic moduli, 
K(q z ), Eq. (6.20), and 

X(q z ) = x(f : y\ (6.29a) 
»{q z ) = ^J\ (6.29b) 

where we have defined the anomalous exponents 

Vn = V\ = -9*F^{x*) , (6-30) 

that take the universal values: 

?7 M = rix ^ 0.6919e , (6.31a) 
= 0.346 for rf = 3, (6.31b) 

controlled by the stable fixed point. The above results 
show that both Lame coefficients /i and A vanish at long 
wavelengths with the same exponent r]^. However, their 
ratio x = \l [i flows to a stable fixed point x* — —0.03272, 
implying a universal negative Poisson ratio 



A 



a P = 



A + 2^ X* + 2 
-0.0166 + 0(e), 



(6.32) 



characterizing the pinned spontaneous vortex solid. 

The above calculations can clearly be extended to other 
directions of q by choosing a more general £* that is a 
function of both q z and q^. To this end we now use a 
more general choice: 

K (t) {q z e e ) A + M (0 (q ± e ur f = n{t)K\ , (6.33) 

determined by the scaling of the denominator of a trans- 
verse part of the phonon correlation function, that for 
q± = reduces to the earlier condition, Eq. (6.17). On 
the other hand, for q z = we obtain the new condition 
(compare to equation (6.17)) 



qj_e 

with the universal anisotropy exponent 
^ = 2 _ + V* 



(6.34) 



(6.35) 



Arguments analogous to the ones developed in the above 
paragraphs show that the above condition, Eq. (6.34), 
leads to the following long-scale behavior of model pa- 
rameters: 



1 

-Vk/C 



(6.36a) 
(6.36b) 
(6.36c) 
(6.36d) 



Combining equations (6.20)-(6.22), (6.29a)-(6.29b) and 
(6.36a)-(6.36d), we finally see that the general scaling 
behavior of the model parameters can be written in the 
form: 

A(q) = x.i*{Qz& L )' , ''U{q±£ L /{Qzg L ) < ) , (6.37a) 
M (q) = n{q x & L )'»>f ll {q±£ L /{q z g L )<) , (6.37b) 
«(q) = <1z^ L )-^Uq^l L /{q z ^ L ) C ) , (6.37c) 
A t (q) = Mq^r^fA^^/iq^^) , (6.37d) 

where we have taken A z = 1/£^ L (which is a reason- 
able choice since is the shortest length scale one 
can go to before perturbation theory can no longer be 
applied) and used the fact that = {^ L f {p/ n) 1 / 2 . 
The crossover functions fn(x), f K {x) and fA t ( x ) have the 
following asymptotic behavior: 



const., 



for 



x v»/C tX -v«/C t x-^t/<, for x » 1. 
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3. Exact relation between exponents 

Since the recursion relation for g, Eq. (6.9a) involves 
functions i^ )K) A t that respectively determined the expo- 
nents rjfj,, k, A t , it quite clearly leads to an exact relation 
between these critical exponents. To see this, we note 
that at the fixed point (dg/d£)\ gt = 0. Combin- 

ing this with the definitions of the critical exponents, Eqs. 
(6.30), (6.24a) and (6.24c), we obtain an exact exponent 
relation: 

(7 - 2d) + f] At = ^ »fo + ^ V. • (6-38) 

From a more fundamental perspective, this exact scal- 
ing relation (Ward identity) is also necessitated by the 
underlying rotational invariance (that is spontaneously 
broken by the vortex lattice). To see this, we recall that 
rotational invariance requires that the form of the non- 
linear strain tensor v a p, Eq. (5.2), be preserved by the 
perturbation theory in these nonlinearities. Performing 
the anisotropic rcscaling in terms of the physical scal- 
ing exponents and imposing this condition (contained in 
<p = 2 — u), we obtain the exact scaling relation, Eq. 
(6.38) above. 

To order e, Eq. (6.38) is given by: 

2e + VA t + \^ = ^ Vk , (6.39) 

and is clearly satisfied by the exponents of equations 
(6.24b), (6.24d) and (6.31b). In d = 3 it reduces to: 

VA t = 2r? K - 1 , (6.40) 

a relation that, in principle, should be experimentally 
testable. 



VII. TOPOLOGICAL STABILITY OF THE 
RAMDOMLY PINNED SPONTANEOUS VORTEX 
LATTICE 

We now turn our attention to the important question 
of whether the new phase of the SV lattice predicted in 
the previous section, and that is controlled by the low- 
temperature fixed point (6.12), is stable against the pro- 
liferation of dislocations. It has been shown that long- 
range translational order in the SV phase is destroyed by 
arbitrarily weak quenched disorder. This does not nec- 
essarily mean, however, that the topology of the phase 
is destroyed. It is the topology, or "connectivity" of the 
phase, that is responsible for the phase's elasticity. If the 
topology is destroyed, then the elasticity is lost. Thus, if 
the SV lattice remains topologically ordered, it will dis- 
play SV solid-like elasticity It will, for example, still have 
a non-zero shear modulus and shear distortions in the x y 
plane will still cost energy. We refer to a phase lacking 
long ranged order, but that still has topological order, as 



an elastic glass. 34,51,52 Of course, the topological order 
is eventually lost if the temperature or disorder strength 
become sufficiently large. For sufficiently large temper- 
ature T or disorder strength A, the topological order of 
the putative elastic glass will be destroyed via a thermo- 
dynamically sharp transition into a translationally and 
topologically disordered phase. 

There are already examples of phase transitions be- 
tween two phases which both lack long-ranged order, the 
most famous being the Kosterlitz-Thouless transition in 
the 2d XY model. 53 In that problem, the transition is 
associated not with the disappearance of long-ranged or- 
der, but rather with the unbinding of neutral pairs of 
topological defects with increasing temperature. The di- 
vergences associated with tilt disorder that destroy the 
long-ranged translational order in the columnar elastic 
glass are of course much stronger (power-law) than in 
the 2d XY model where they are logarithmic. There are, 
however, examples of phases in much more strongly dis- 
ordered systems in which topological defects nonetheless 
remain bound. 54 The possibility of such a transition from 
a translationally disordered smectic-A liquid crystal has 
also been considered by Radzihovsky and Toner. 29 

In general, the destruction of topological order in a 
phase occurs via the proliferation of defects. In the 2d 
XY model these defects are spin vortices. In the SV 
lattice the main mechanism for the destruction of topo- 
logical order is through the proliferation of dislocations. 
An edge dislocation, which is an insertion of an extra 
row of flux lines is shown, for a square lattice, in Fig. 5. 
The Burgers vector (the extra lattice vector) associated 
with this dislocation is b = ay. In the remainder of this 
Section, we restrict ourselves to spatial dimensions d = 3. 

This edge dislocation will reduce the elasticity against 
shears in the yz plane. One could also have an edge dislo- 
cation corresponding to multiple extra inserted rows and 
in this case the Burgers vector will be a multiple of the 
lattice vector. If the rows of vortex lines are inserted in 
the yz plane then the Burgers vector will be an integer 
number of lattice vectors in the x direction. One could 




Dislocation 



FIG. 5. A square spontaneous vortex with an extra row of 
flux lines introduced from the right, along the x axis. The 
core of the dislocation comes out of the page. 
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also have both rows inserted in the yz and the xz planes, 
and the corresponding Burgers vector would be a sum of 
lattice vectors. Of course, one can also consider dislo- 
cations in a hexagonal lattice, the only difference being 
that the lattice vectors do not simply point along x and 

y- 

When the displacement field u(r) is integrated around 
a closed curve enclosing the core of the edge dislocation 
shown in Fig. 5, the result is not zero as it would be 
without a dislocation, but rather b = a 2 = ay, i.e., 

Idu = a 2 . (7.1) 

In general, a dislocation is a line or curve, with the prop- 
erty that an integral of u around any curve enclosing the 
dislocation core is some Burgers vector b = Mai + Na 2 , 
i.e., 

j>dvL = h = Mai + Na 2 , (7.2) 

where M and N are integers. In the cases of square 
and hexagonal lattices, in which the lattice vectors are of 
equal length, Eq. (7.2) can be rewritten as: 

<j> dui = Ma , (7.3a) 

<j> du 2 = Na , (7.3b) 

where Ui is the i'th component of u, defined in the co- 
ordinate system with (not necessarily orthogonal) unit 
vectors ai = ai/a, k 2 — & 2 /a, i.e., u = U\ai + u 2 a 2 . For 
a square lattice i is either x or y. 

Eqs. (7.3) apply for all types of dislocations. The edge 
dislocations, discussed above, have cores that run parallel 
to the vortex lines but in general dislocation lines can run 
in any direction. An extreme case is a screw dislocation 
which runs perpendicular to the vortices. 



A. The Defect Elastic Hamiltonian 

Eqs. (7.3a)-(7.3b) imply that in the presence of a dislo- 
cation the u field is no longer single valued. Given that a 
single valued u field is obviously a prerequisite of a mean- 
ingful elastic theory it is necessary to introduce a "cut" 
in the following fashion. In order to remain single valued 
as the core of the dislocation is circulated, there must be 
points at which the value of u jumps discontinuously by 
an amount equal to the Burgers vector, b. The cut is the 
locus of the points at which this jump occurs, and the 
value of u on either side of the cut differs by b. In the 
absence of dislocations, the elasticity of the system can 
be described using a field 

v 4 = Vu, , (7.4) 



which is continuous and smooth almost everywhere. 
However, the presence of a cut means that Vj will be 
singular at the cut. The cut, however, does not represent 
anything physical, and therefore, nor does a singular Vj. 
There is nothing singularly elastic happening in the real 
system as the (arbitrarily chosen) cut is crossed. There- 
fore, a modified Vj, from which the unphysical, singular 
piece is removed, is used to construct an elastic theory 
in the presence of a dislocation. The modified is equal 
to V«i everywhere except at the cut. 

Using Stokes law, Eqs. (7.3a)-(7.3b) can be rewritten 
in differential form: 

V x v, = m 4 . (7.5) 

Since v, is no longer the gradient of a single valued field, 
it can have a non-zero curl. This is where the presence 
of the dislocation is manifest. In the absence of a dis- 
location the curl would be zero. The field m is defined 
as: 

m 4 (r)=^ f aiNjtjir^s^ir-Ti^dsj , (7.6) 

3 

where sj parameterizes the j'th dislocation loop, Vj(sj) 
is the position of that loop, tj(rj(sj) is its local unit tan- 
gent, and Nj the "charge" or number of excess lattice 
vectors along a^ associated with the dislocation. Note 
that Nj is independent of Sj , since the charge of a given 
line is constant along the line defect. It is important to 
note that Vj is not the i'th component of the three dimen- 
sional vector v. Rather, it is one of the two vectors vi 
and v 2 , constructed using Eq. (7.4). The same applies 
to mj. 

Eq. (7.5) implies that 

V -111,(0=0, (7.7) 

which simply means that dislocation lines cannot end in 
the bulk of the sample: they must either form closed 
loops or extend entirely through the system. It is useful 
to work with Fourier transformed fields: 

v»(q) = iq«»(q) , (7.8a) 
m i (q)=iqxv,(q), (7.8b) 

and work with the longitudinal and transverse compo- 
nents vi,(q) and Vr(q), defined by: 

vi,(q) =iqu L (q), (7.9a) 
v T (q) = iqu T (q) , (7.9b) 

where Uz,(q) and iir(q) are the longitudinal and trans- 
verse components of u(q) given, in d = 3, by 

u(q) = it L (q)qu + u T (q)(z x qj_) . (7.10) 

The field Vj can be easily constructed from v^(q) and 
vr(q) according to: 
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Vi(q) = v L (q)(q±)i + u T (q)(z x q±) t 



(7.11) 



where (q_i_)i is the z'th component of q_i_, and (z x q_i_)j 
is the i'th component of i x q^. Similarly, 



mj(q) = mi,(q)(qx)i + m T (q)(z x qj.)* 



(7.12) 



In Fourier space the condition that dislocation lines must 
form closed loops or extend entirely through the sample, 
Eq. (7.7) gives: 



q- m L (q) =0, 
q m T (q) = 0, 



(7.13a) 
(7.13b) 



and the boundary condition, Eq. (7.5) corresponding to 
a topological defect is: 



iq x vz, (q) = m L (q) , 
iq x v T (q) = m T (q) . 



(7.14a) 
(7.14b) 



The real space Hamiltonian to harmonic order in d = 3 
dimensions is given by Eq. (5.3) that we rewrite here for 
dcfinitcncss: 



H[u(r)} 



dr — 

2 



h(r) • d z u] , (7.15) 



where the reader is reminded that h(r) is the random tilt 
disorder, with short-ranged, isotropic correlations 



Mr)MrO = A t <5 a/ j5(r - r 7 ) . 



(7.16) 



The above real space Hamiltonian can then be rewrit- 
ten in Fourier space, where u(q) and h(q) can both be 
decomposed into longitudinal and transverse pieces ac- 
cording to Eq. (7.10). The longitudinal and transverse 
modes decouple and one obtains 

tf[u(q),h(q)] =H L [u L (q),h L (q)]+H T [u T (q),h T (q)] , 

(7.17) 



with 

H l/ t[u l /t{Ql), h L/T (q)} = \Y1 [ B l/t (q\ + ^l/rlt) K/r(q)| 2 - iq z u{q)h L/T {-q) 



(7.18) 



where Bl = A + 2/z , Bt = /U and \\j T = n/B L / T . This 
Hamiltonian does not have any anomalous elasticity in 
it. It is instructive to perform the dislocation theory at 
the harmonic level first and then to generalize to a full 
anharmonic theory afterwards. 

Using the decoupled Hamiltonians for Ui(q) and Ur{q) 
it is possible, using the constraints Eqs. (7.13a)-(7.13b) 
and (7. 14a)- (7. 14b) to obtain decoupled elastic Hamilto- 
nians for the dislocation fields mj;(q) and m-r(q). The 
first step is to minimize it [u l it (q) j^l/t (q)] with 
respect to u L / T (q). The Euler-Lagrange equation for 
Wl/t(q) is given by: 



( ? i + A| /T ^). vr( q)-^M=0, (7.19) 



and can be rewritten in terms of /t (q) as 

q z h L /T(q) 



q.L • v L/T (q) + \ L/T q z v z L/T {q) + 



B 



L/T 



(7.20) 



where , T (q) and v z L ^ T (q) are the _L and z components, 
respectively, of v^/y(q). The next step is to solve Eq. 
(7.20) above and Eqs. (7.14a)-(7.14b) for v i/T (q). Do- 
ing this, one obtains: 



vl/t(q) 



iq x m L/T (q) q z q f i(m L/T (q)) 



-\L/T 



-(i-a! /t ^ 2 ) 



frz,/T(q) 



B 



L/T 



(7.21) 



where T^ T is defined as 



q — l_L A L/Tl. 

These expressions for v^ T (q) can then be inserted 



which gives the defect Hamiltonian: 
(7-22) „._„Li„_ f ^ ,.,„m , uTf 



H d = H^[m L (q),h L (q)}+Hi [m r (q), Mq)] , (7.23) 



back into the original elastic Hamiltonian of Eq. (7.17), where 



H d /T [m L/T (q), h L/T (q)} = J2 



Kq^\m 



L/r(q) 



B 



L/T\m L/T 



(q) 



^L/T 



m L/T (q) • a L/T (-q) 



(7.24) 
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and where & L / T (q) is related to h L / T (q) via: 
ih L /T(ci)(z x q) 



ai/T(q) = 



(7.25) 



B. Duality Transformation 

The defect Hamiltonian of Eq. (7.24) makes it possible 
for us to study the statistical mechanics of the disloca- 
tions. We do this by performing a duality transformation. 
As will be seen below, to perform this transformation it is 
necessary to have a partition function in which the trace 
is taken over the real space dislocation fields. However, 
as will also be seen below, although this is possible, it 
is not necessary to rewrite the defect Hamiltonian above 
(which is expressed in terms of the fields m L / T (q)) in 
terms of the fields nil (q) and m 2 (q) . The reader is re- 
minded that mi(r) and ni2(r) are the dislocation loop 
fields associated with dislocations along ai and a 2 respec- 
tively, where ai and a 2 are basis vectors of the lattice. 

Henceforth for simplicity a square lattice will be as- 
sumed, i.e., ai = ax and a 2 = ay. This assumption 
makes certain parts of the calculation more straightfor- 
ward than they would be for arbitrary ai and a 2 , e.g., for 
a hexagonal lattice. This assumption may seem invalid 
given that the Hamiltonian that has been used thus far is 
that of a hexagonal lattice. However, the determination 
of topological stability is ultimately based upon scaling 
arguments. Since the scaling, e.g. of the propagators 
and correlation functions, is the same in the hexagonal 
and square lattices the answer, which depends only on 
this scaling, obtained using this innocuous simplification 
to a square lattice should be the same as that for the 
hexagonal lattice. This assumption would not, however, 
be so innocuous for a theory of the transition, where the 
structure of the lattice is important. 

The first step in the duality transformation is to put 
the model on a simple cubic lattice (to make it well- 
defined at short distances); now, nii(r) and m 2 (r) are 
defined on the sites r of the lattice, and they take on 
values 



m i (r) = ^(n?(r),n?(r),n?(r)) ) 



(7.26) 



where nf is the a'th component of and the Einstein 
summation convention has been suspended. Each of nf, 
n\, and n\ is an integer. For the square lattice considered 
here the lengths of the two basis vectors are the same so 
that a\ = a 2 = a. d is the cubic lattice constant used in 
the discretization. The partition function for this model 
is then 



Z[{h}} = ]T 



E < 

Hi(r)} 



-S[{mi}] 



(7.27) 



5[{mJ]=/3[i/ d [{m 1 ,m 2 }]+i? c ^|m 4 (r)| 2 j , (7.28) 



and the sum is over all integer-valued configurations of 
nij's given by Eq. (7.26), satisfying the dislocation line 
continuity constraint (which is denoted by the prime): 



V • rrij = 



(7.29) 



where the divergence now represents a lattice divergence, 
Hd is obtainable from Eq. (7.24), and a core energy term 
£ c ^ r i |nii(r)| 2 has been added, to account for energies 
near the cores of the defect line that are not accurately 
treated by the above continuum elastic theory. The sym- 
metry of the square lattice dictates that the core energies 
be the same for both rrii(r) and m 2 (r). The reader's at- 
tention is directed to the fact that the partition function 
still depends implicitly on the configuration of the ran- 
dom fields {h}. 

To proceed, the constraints V • nil = and V • m 2 = 
are enforced by introducing new auxiliary fields Q\ (r) and 
2 (r), and rewriting the partition function Eq. (7.27) as: 



z = IiJ de 1 (v)de 2 (v) 



E E < 

i {m;(r)} 



-S[W.}]+^ r , e s (r)V-m t (r)<i 2 /a 



(7.30) 



where the sum over {mi} is now unconstrained (and no 
longer has a prime), the constraint V- = is enforced 
by integration over 9i, since 



S(V -111,)= / <Wi(r)e 
Jo 



ie,(r)V-mi(r)d 2 /a 



(7.31) 



where the 5 is a Kronecker delta since m^, and, hence, 
V • nij, is integer valued. 

Now one can "integrate" (actually sum) by parts, and 
rewrite 

]T 0i(r)V -mdr) = -J2 m 4 (r) • W;(r) 

r,i r,i 

+ surface terms . (7.32) 

The next step is to introduce dummy gauge fields Ai and 
A 2 to mediate the long-ranged interaction between the 
defects loops mi and m 2 in the Hamiltonian Eq. (7.23). 
This is accomplished using the Hubbard-Stratonovich 
transformation and rewriting the partition function as: 



Z = ]J d6 1 (r)de 2 {r)dA 1 (r)dA 2 (r) 

r •> 

X E E e~ s ^ m ^ e ^5{V • A,) 

i {m,(r)} 



(7.33) 



where 



with 
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5 = k E { mi(r) • [ - ^ V ^( r ) + rf3 (^( r ) + a *( r ))] + ^c|m 4 | 2 } 

q 



(7.34) 



where A^(q) and Ay(q) are defined in the same way as 
vl/t(q) and m L/T (q), i.e., 

A L (q) - qiAi(q) + q 2 A 2 (q) , (7.35a) 
A T (q) = (z x q)iAi(q) + (z x q) 2 Ai(q) , (7.35b) 

where (q_i_), is the i'th component of q^ and (z x q_ L ) i 
is the i'th component of z x qi, and where a,(r) is the 
Fourier transform of 



a»(q) = a L (q)(q ± ) 4 + a T (q)(z x qu), . 



(7.36) 



The two goals of all of these manipulations have now 
been achieved: the sum on {rrii(r)} is now unconstrained, 
and the sum on each site over mj(r) is now decoupled 
from that on every other site. Furthermore, this sum is 
readily recognized to be nothing more than the "periodic 
Gaussian" first used by Villain. 38 The partition function 
Eq. (7.33) can thus be rewritten: 



= 11 / <»i(r)^2(r)dAi(r)dA 2 (r)5(V-Ai(r))<5(V-Ai(r)) 
x exp J - W( r + *a) - W - A ?( r ) + <(')] 

I, r,i.a 



:|Af(q)| 2 + -f |Al(q 



Kq 



<L| A ^(q)| 2 + ^|A^' 2 



g 



(7.37) 



where Af (r) is the a'th component of A, (r) and the well- 
known 27r-periodic Villain potential V p (x) 



-v P (x) 



E 



2 E c /T+ixn 



(7.38) 



has the usual property that the smaller E c /T is (i.e., the 
higher the temperature in the original random-tilt smec- 
tic model), the sharper the potential minima. This can 
be seen by looking at two extremes. When E c /T is zero, 
the Villain potential is a just a periodic set of delta func- 
tions. When it is very large, then the n = 0, ±1 terms of 
the series dominate and the potential is a constant plus 
a cosine of x which is obviously not as sharp as a set of 
delta functions. Furthermore, the amplitude of the co- 
sine is proportional to e~ E "l T and hence gets smaller as 



t is reduced. Thus raising the temperature in the original 
model is like lowering the temperature in the dual model 
Eq. (7.37). It is precisely this familiar temperature in- 
version associated with duality that leads to an inverted 
XY transition for bulk smectics 55 and three-dimensional 
disorder- free charged supcrfluids 56 . It also plays an im- 
portant role here, as we shall see in a moment. 

Standard universality arguments imply that replacing 
the periodic potential V p {x) in Eq. (7.37) by any other 
non-singular periodic function should not change the uni- 
versality class of the transition. In particular, we could 
replace V p (x) by cos(x). The resultant model would be 
precisely the "fixed length" version of the "soft spin" , or 
Landau-Ginsburg- Wilson model, with the complex "ac- 
tion" 



{r} 



ad 



+ 



\ E l A L(q)l 2 + ^|Ai(q)| 2 + 5l A T(q)l 2 + 5lA^(q)| 2 



2 '■ — ' inq 
q 



Bt 



(7.39) 



where Vi( r ) is a complex field whose phase is 0j(r), and 
the reduced temperature t and quartic coupling u are pa- 
rameters of the model (and are the same for tpi(r) an d 
^(r)) because of the symmetry of the square lattice). 
Because of the duality transformation's inversion of the 
temperature axis, the reduced temperature t is a mono- 



tonically decreasing function of the temperature T (of 
the original dislocation loop model), which vanishes at 
the mean- field transition temperature Tmf of the fixed 
length model Eq. (7.37). 

Universality also implies that this "soft-spin" model 
should be in the same universality class as the fixed 
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length model Eq. (7.37). Therefore model Eq. (7.39), 
shall henceforth be used because it is more straightfor- 
ward to analyze perturbatively. 

Before undertaking that analysis, it is worth noting 
another consequence of duality inversion of the tempera- 
ture axis: the ordered phase of the dual model Eq. (7.39) 
corresponds to the disordered (i.e., dislocation loops un- 
bound) phase of the original dislocation loop gas model. 
That is, the low dual-temperature phase described by 



Wr)> + 



(7.40) 



corresponds to the topologically disordered, dislocation 
unbound phase of the SV lattice. 



C. Analysis Of The Dual Model 

Disorder is included in Eq. (7.39) through the 
quenched gauge-fields ai(r) and a 2 (r), which are related 



to the random tilt field h(r) by Eqs. (7.25)-(7.36). The 
partition function 

Z[{h}} = J [# 1 # 2 ][rfA 1 rfA 2 ]e- s ^' A *- h ] 

x <J(V • Ai)J(V • A 2 ) , (7.41) 

with S given by Eq. (7.39), is thus an implicit function 
of the random tilt field configuration {h(r)}. 

The dependence of the partition function Z on the 
quenched field {h} is dealt with by using the replica trick, 
Eq. (5.9). Doing this leads to the replicated partition 
function: 

/n 
[dsn] [da 2 ] J] [#i«] [dAi„] [dA 2a ] 

o=l 

x e- s ^' A ^P[ ai }P[a 2 }S(V • Ai)5(V ■ A 2 ) , 

(7.42) 

with: 



r a— 1 



faa 



^(iA iQ +a 4 ) fa* a + t\fa a \ 2 + u\fa, 



| 4 + «#ia| 2 hM 2 } 



+ 9 E b^l A i«(q)l 2 + ^IAL(q)! 2 + ^Ta(<D\ 2 + ±\ A Ta(<l)\ 



Bt 



(7.43) 



Note that A la is the a'th replica of the A x field, not the ath component of Ai. The probability distribution P[aj] of 
the field in Eq. (7.42) is Gaussian, defined by Eq. (7.16). Thus the distributions of P[ai] and P[a 2 ] are completely 

specified by the averages (4(q)aj(— q) and a 2 (q}a 3 2 (— q) respectively. These are easily evaluated using the relations 
(7.25) and (7.36). For ai: 



a i(qK(-q) = ai(q) a i(-q)((q-L)i) + afr{q)a J T (-q)((z x q±)i) 2 + 2a J L (q)a^(-q)(q±)i(z x qj_)i 



((g^r) 2 ((Sxqi)i) 

(r£) 2 



(IT) 2 



(7.44) 



where the second equality was obtained using Eq. (7.16), 

SizSjz- Similarly, 



and where we defined 5^ = S t] 



4(qK 2 (-q) = <^A t<? i 



((q±) 2 ) 2 ((zxq^F 



(7.45) 



One now must consider the statistical mechanics of the 
model defined by Eqs. (7.42) through (7.45), in the limit 
n -> 0. 

Since it is only necessary to determine whether the 
topologically ordered columnar elastic glass phase is sta- 
ble against disorder induced dislocation loop unbinding, 
a complete analysis of the critical properties of the dis- 
location loop unbinding is unnecessary. Thus, it is suffi- 
cient to follow the example of Radzihovsky and Toner, 29 
and consider only the limit in which As a result, 

the fa fluctuations are subdominant to to those of the 



gauge field Ai allowing a mean field treatment of fa. In 
this case the remaining gauge fluctuations may be treated 
exactly. The stability of the columnar elastic glass is as- 
sessed by studying the effect of diagrammatic corrections 
on the reduced dual temperature t. The lowest order con- 
tribution to the renormalizcd dual temperature tn comes 
from the average of the "diamagnetic" terms 




+ 



FIG. 6. The Feynman diagrams that dominate the renor- 
malization of the reduced dual temperature t. 
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^=^EE(<N 2 )-W 2 )W 2 , (7.46) 



graphically illustrated in Fig. 6. 

Taking into account the symmetry of the square lat- 
tice, and first angular averaging in the _L plane, the cor- 
rection to the renormalized reduced dual temperature is, 
in d = 3 dimensions, 



tR — t() + 



ca 2 d 2 

J^2 



K-q z 



1 

g 



1 

q 



Bl 

T L 

q 



-A t gi 



9 



(7.47) 



where the first four terms proportional to l/Tq^ T come 
from the first graph in Fig. 6, with the wiggly line repre- 
senting Aj fluctuations, while the terms proportional to 
At come from the second graph with the internal dotted 
line representing the quenched gauge field a^. The nega- 
tive sign of the disorder contribution leads to an increase 
in the dual T c and can be traced back to the fact that 
the action S r in Eq. (7.43) is complex-valued. 

Performing this integral, one finds no infrared diver- 
gence in spatial dimensions d = 3 (the only dimensions 
in which the analysis of this section applies), indicating 
that the renormalized t r remains finite as the system size 
diverges. A finite t r implies that there is a temperature 
regime in which the dual parameter is in its disordered 
phase, which in turn implies that there is a temperature 
regime in which dislocations in the columnar phase re- 
main bound. Thus, in the harmonic model, there is a 
stable columnar elastic glass. 



fields mj(r), and hence, those fields would not be decou- 
pled by a simple Hubbard-Stratonovich transformation, 
as they can in the harmonic model. 

For all these reasons, a completely honest treatment of 
dislocations in the full, anharmonic model is extremely 
difficult. The next best, tractable, (but uncontrolled) 
approach, introduced by Radzihovsky and Toner, 29 is to 
simply replace the elastic constants /it, A and k, and the 
tilt disorder variance A t in Eq. (7.47) for the renormal- 
ized dislocation unbinding transition temperature, with 
the renormalized wavevector dependent moduli yu(q), 
A(q), n(q), and At(q), derived in Sec. VI B. One can 
hope to justify this procedure by doing an RG matching 
calculation. 

Doing this, simple power counting shows that in the 
IR limit (i.e., at long length scales) it is the negative, 
disorder, contribution to the renormaliztion of the dual 
temperature which dominates, i.e., 



St = -0(1) x 



A t (q)«i 



T 2 ,/ q (<zi + A!(q)g z 4 ) 2 



(7.48) 



where A t (q) is given by Eq. (6.22) and A^(q) = 
K(q)//i(q). The contribution of the piece proportional 
to l/(r^) 2 is absorbed into the 0(1) coefficient, since 
in the IR limit and scale the same way with q 
(because /Lt(q)/A(q) is independent of q). Imposing an 
infrared cutoff q z > L _1 on the wavevector integral in 
Eq. (7.48) one finds that the infrared piece scales with 
the linear size of the system L like 



St 



IR 



L V ca 2 d 2 A t 



(7.49) 



D. Treating The Nonlinearities 

Of course, one would like to also perform a full anhar- 
monic theory of dislocations (since the anharmonic ef- 
fects are obviously important) but unfortunately such a 
theory is simply intractable. In particular, the fact that, 
in an anharmonic theory, the interaction energy between 
dislocations cannot be written as a sum of pairwise in- 
teractions (since their fields do not simply add) makes it 
impossible to even write down a general expression for 
the energy of an arbitrary dislocation configuration. At 
best, one might hope to be able to write down the energy 
for a few simple, high symmetry configurations (e.g., a 
simple, straight dislocation line). Such specialized results 
would be of no use in a full statistical theory of defect 
unbinding, which requires consideration of very compli- 
cated, tangled configurations of dislocations, which for 
cntropic reasons dominate the free energy near the dislo- 
cation unbinding transition. 

Furthermore, even if one could write down the anhar- 
monic energy for an arbitrary dislocation configuration, 
it would presumably be anharmonic in the dislocation 



with 

7 = VA - 1 ■ (7.50) 

If 7 > the negative correction to t diverges, correspond- 
ing to an infinite downward renormalization of the dual 
temperature, even for arbitrarily weak disorder, implying 
that dislocations in the SV lattice are always unbound. 
On the other hand, if 7 < then there is no infinite neg- 
ative divergence of the dual temperature implying that 
the columnar elastic glass is stable in a finite range of 
temperatures. 

Using the exact scaling relation in d — 3, Eq. (6.38) 
allows 7 to be reexpressed in terms of the more experi- 
mentally accessible n K exponent. Doing this one finds 

7 = 2( Vk - 1) . (7.51) 

Thus for n K < 1 the columnar elastic glass phase is sta- 
ble. In d = 3, the estimate of n K , found in the earlier 
e = 7/2 — d expansion, is n K = 0.74, so that a colum- 
nar elastic glass is indeed predicted to be stable against 
dislocation unbinding. 
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FIG. 7. The region of w K and r? M that satisfy the criteria 
Eqs. (7.52a)-(7.52c) for a long-ranged orientationally ordered, 
stable (at sufficiently low disorder strength) "columnar elastic 
glass" phase. 

Thus, in summary, for a topologically ordered phase to 
exist the anomalous exponents must satisfy three distinct 
criteria, corresponding to a destruction of long-ranged 
translational order, the existence of long-range orienta- 
tional order and stability against dislocation unbinding. 
These three criteria are respectively: 



»?«<!, 



(7.52a) 
(7.52b) 
(7.52c) 



and are along with a fourth bound that r\ K > 1/2 (which 
is equivalent, via the exact scaling relation in d — 3, to 
the condition that r)A t > 0), summarized graphically, in 
Fig. 7. The exponents found in the e = 7/2 — d ex- 
pansion satisfy all of these bounds; hence a topologically 
ordered phase of the SV lattice is possible in three di- 
mensions. We shall term this phase (which is analogous 
to the Bragg glass phase of vortices in disordered type 
II superconductors 34 ) a "columnar elastic glass" (CEG) 
phase. 



VIII. MAGNETIZATION CURVE AND OTHER 
PHENOMENOLOGY 



We expect that for length scales longer than £ NL the 
strong, power-law anomalous elasticity will manifest it- 
self in the behavior of all physical observables of the SV 
solid, such as, e.g., the width of the structure function 
peak 5*(q) measured in neutron scattering and the be- 
havior in its tails. Anomalous elasticity also implies that 
the stress (cr)-strain (e c ) relation becomes non-Hookean, 
i.e., nonlinear, even for arbitrarily small applied stress 
a. (Note here, that a and e c should not be confused 
with the Poisson ratio, <r p and e-expansion parameter, re- 
spectively). In this section wc will discuss the expected 
behavior of the magnetization curve and the structure 
function. 



A. Nonlinear Magnetization Curve 

The origin of the non-Hookean stress-strain behavior 
can be seen most clearly for a purely compressive stress, 
dij = a5ij , (which, as we shall see shortly, arises from an 
externally applied magnetic field H along the flux lines). 
Such a stress adds a term proportional to 



-crV • u, 
and a term proportional to 
<r{d z n) 2 



(8.1) 



(8.2) 



(which explicitly breaks rotational invariance) to the 
original Hamiltonian. The actual values of the dimen- 
sionless constants of proportionality are not important 
and for simplicity we set them equal to unity. In Fourier 
space, the second, symmetry breaking term becomes 

^ 2 |u(q)| 2 , (8.3) 

and begins to dominate over the K(q)g 4 |u(q)| 2 term once 

aq\ > K (q)g 4 . (8.4) 

This clearly happens for q z 's smaller than a critical q c 
given by 



c<?c = K (Qz = q c ,q± = 0)q 4 c 



.5) 



Using equation (6.20) for /c(q) in this expression and solv- 
ing for q c , we find 



where we denned 

5 = 2 - Vk . 



(8.6) 



•7) 



The important thing to note about this result is that it 
(singularly) depends on the applied stress a. 

Now, for q z > q c (a) the new, stress- induced rotational 
symmetry breaking term crg5;|u(q)| 2 in the Hamiltonian is 
negligible (subdominant) relative to the vortex curvature 
energy K(q)q 4 |u(q)| 2 . However, on longer length scales, 
qz <C 5c(f), the er-term dominates, suppressing the SV 
solid fluctuations and hence cutting off the anomalous 
elasticity. Therefore, the elastic moduli saturate, for 
qz -C q c (c) , at their values at q z — q c {cr), q± = 0, given 

by 



A(q 



0) 
0) 



Klz 

Kiz 



g c ,q_L 



o), 
o), 

Jill 



(8.8a) 
(8.8b) 



and, hence, both are proportional to qT ■ Since q c is a 
function of the applied stress a, as shown in equation 
(8.6), this implies that the long wavelength limit of the 
2d bulk modulus: 
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C(a) = fj,(q c ) + X(q c ), 
™ n[qc(<r)g L ] r "' , 

-)'■ 



with a < ctatl = L )- 2 +^-^/^, and 



/3 = 



2-r?« 



(8.9) 



(8.10) 



To obtain the response of the system to a uniform, 
isotropic compressive stress within the plane transverse 
to the spontaneous vortex axis, the 2d bulk modulus C 
times the compressive strain e c = V • u must be balanced 
against the applied compressive stress, i.e., C(cr)e c = a. 
Hence, we find that for a stress a < <jnl, the SV glass 
exhibits nonlinear (non-Hookean) strain-stress relation: 



e c (a) = a/C{a) . 
occr", 



with 



a 



1-/3, 

= 1 - 
2 

w 0.72, 



(8.11) 



(8.12a) 
(8.12b) 



where in the last equality, the numerical 3d values of rj^ 
and r\ K obtained from the e-expansion have been used. 

To see that the non-Hookean stress-strain relation, Eq. 
(8.11), leads to a nonlinear relation between B and H, 
note that the presence of an external field, H, necessi- 
tates the inclusion of an additional B H term in the 
Hamiltonian. A field H along the z direction (the direc- 
tion of the initial magnetic induction B(iJ = 0)), will 
therefore induce an increase in the magnetic induction, 
SB(H) = B(H) - B(H = 0). In the vortex lattice state, 
however, since each flux line carries exactly one flux quan- 
tum, (f)Q, the only way to increase B is to move the vortex 
lines closer together, i.e., compress the lattice. This im- 
plies that the external field leads to a compressive strain 
term, proportional to — if V- u. The external field H also 
explicitly breaks rotational symmetry, leading to the ad- 
dition of a term ~ H\d z u\ 2 . Thus, the application of a 
field H leads to a compressional strain (proportional to 
5B) and to explicit rotational symmetry breaking, just 
as a uniform compressive stress does. The two cases are 
exactly analogous, with a replaced by 57 H, and the non- 
Hookean stress-strain relation, Eq. (8.11), immediately 
translates into a nonlinear SB(H) — H relation of exactly 
the same form; i.e., into 



5B(H) cx H a , 



.13) 



with the universal exponent a = 0.72 ± 0.04. This rela- 
tion is shown in Fig. 1. 



Now, of course, in a real crystalline ferromagnet, 
there are always crystalline symmetry breaking fields, 
which pick out a preferred direction for the magneti- 
zation, thereby explicitly breaking the rotational invari- 
ance. These symmetry breaking fields lead to a non-zero 
tilt modulus, i.e., to a term 



-V x \d z uf 



.14) 



in the Hamiltonian, which cuts off the anomalous elas- 
ticity even in the absence of an external field H. Con- 
sequently, we expect the predictions of our rotationally 
invariant theory to break down beyond a crossover length 
scale 



(f ) 1/S (^ L ) 



-Vk/S 



(8.15) 



If the crystal symmetry breaking fields are weak, V x will 
be small, and our anomalous elastic theory will be valid 
for a wide range of length scales L z ,± satisfying 



(8.16) 



For longer length scales the soft anisotropic elasticity of 
the SV lattice will cross over to the conventional "ten- 
sion" elasticity of ordinary vortex lattices. This crossover 
in length scales manifests itself in a crossover value H x 
for the applied magnetic field H below which the non- 
Hookean law, Eq. (1.2), breaks down, to be replaced by 
a conventional, Hookean relation SB cx H, as illustrated 
in Fig. 1. The value of the crossover field H x is propor- 
tional to the crystal symmetry breaking field V x ■ 

The order of magnitude of H x is determined by V x , 
which in turn can be estimated from the width of the 
hysteresis loops for the magnetic superconductor in the 
spontaneous vortex lattice state as V x ~ M H , where 
M is the remnant (spontaneous) magnetization and H 
is the coercive field. While we were not able to find any 
hysteresis data of ferromagnetic superconductors in the 
literature, we hope that our work will stimulate such ex- 
periments as well as other attempts to estimate H x , such 
as microscopic models and numerical methods. 

The presence of the crystal symmetry breaking fields 
also makes it experimentally possible to apply the ex- 
ternal magnetic field opposite to the direction of the 
spontaneous magnetization. While the true ground state 
clearly has a "spontaneous " magnetization which is par- 
allel to the externally applied field, this anti-parallel con- 
figuration is naively metastable if |H| < V x . Writing 
H = —Hz, we readily see that, when H reaches the 
metastability limit V x , the iV><|<9 z u| 2 term is exactly 
canceled by the — ^H\d z u\ 2 term arising from the H • M 
interaction. Thus, at precisely this metastable limit, the 
flux lattice becomes "soft" in almost exactly the sense 
we have been discussing throughout this paper. There- 
fore, we expect anomalous elasticity out to arbitrarily 
long length scales at this metastable point, even in the 
presence of crystal symmetry breaking fields. 
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Of course, the exact metastability limit can never be 
reached; as it is approached, the metastable barrier gets 
smaller and smaller, and before it can be reduced to zero 
the system will thermally tunnel over the barrier. Thus, 
it might appear that the point of exact "softness" can 
never be reached, and, hence, that there is no regime 
with anomalous elasticity out to arbitrarily long length 
scales. 

However, this is not necessarily the case. Recently, it 
has been found 32 that if the crystal symmetry breaking 
field has the right mixture of higher spherical harmonics 
in its dependence on the orientation of the flux lattice, it 
is possible for the metastable state to remain metastable 
right up to the point of perfect "softness ". In a sense, 
the metastability actually persists beyond this point: in 
fact, as we move beyond the point of perfect softness, 
the coefficient (call it D(H)) of |<9 z ii| 2 in the Hamilto- 
nian becomes negative, inducing an instability not to the 
true ground state with the flux lattice tipped 180° rela- 
tive to the metastable state, but, rather, to a local mini- 
mum very near the original metastable state which moves 
continuously away from the old metastable position as we 
continue to make D(H) more negative. Indeed, this looks 
very much like the tilting of the molecules away from the 
layer normals in a smcctic A to smectic C transition. We 
will discuss this fascinating and novel transition more in 
the Sec. IX below. 

It should also be noted that many of the ferromagnetic 
superconductor materials are characterized by strong 
easy-plane anisotropy and only much weaker in-plane 
crystal-symmetry breaking fields. Our theory can eas- 
ily accommodate such systems by confining fluctuations 
of the spontaneous vortex lattice to lie in a plane (that we 
define to be the xz plane) by setting u y — 0. The symme- 
try breaking field V x then describes the remnant in-plane 
(xz) anisotropy. The remnant anisotropy and associated 
V x should be quite small in such layered easy plane ma- 
terials. In this case the system is in the universality class 
of the "hybrid columnar Bragg glass" phase described in 
Ref. 58 , which mixes the properties of an "to = 1 smectic 
clastic glass" 31 with those of a conventional Bragg glass. 



B. The Spontaneous Vortex Lattice Structure 
Function 

Even though the neutron scattering peaks of the SV 
lattice will be broad (due to the absence of long-ranged 
translational order) , they will still contain a great deal of 
information about the anomalous elasticity of the CEG 
phase. Since, as demonstrated in Appendix VII (see also 
ref. 21) the CEG phase possesses long-ranged orienta- 
tional order, single domain, rather then powdered aver- 
aged, neutron scattering should be possible. 59 Thus, the 
anisotropic scaling information which is usually lost in a 
powder averaged x-ray scattering experiment would be 
retained, allowing detailed tests of the quantitative pre- 
dictions for r] K , i]^ and rj\ t . 



The neutron scattering will have broadened spots in 
the _L plane centered at the reciprocal lattice vectors. 
The intensity of a Bragg spot centered at reciprocal lat- 
tice vector Q is: 



/(q) = (|PQ(q)l) 2 , 



(8.17) 



where pq (q) is the Fourier transform of piece of the spa- 
tially varying density of vortex lines modulated along Q. 
For any lattice u can be defined via 



PQ(r) = |PQ|e jQ ' (r+u(r)) 



and so the scattering intensity of Eq. (8.17) can be ex- 
pressed as: 



J(q) = | / g Q | 2 (^| el Q-(u(r)-u(0)) e4 (Q-q). 



\P Q \ 2 I e *(Q-i)-' e -<(Q-(»W-«e») a >/2 , (8.19) 



where in going from the first to the second line the 
property of a Gaussian distributed variable x, (e tx ) = 
e -( x >/ 2 ; h as been used. The effect of the u dependent 
exponential is to broaden the peak at Q. This exponen- 
tial can be expressed in terms of the correlation function 
C a p(v) = (Mr) - u a (O)][u (r) - 1^(0)]): 



e -<(Q-(u(r)-u(0))2)/2 = e -Q a Q0C a p(r)/2 



(8.20) 



From Eq. (8.20) the widths in the q z and q± directions of 
the broadened peaks can be obtained. The peak widths 
in the q z and q± directions are defined as the inverse of 
the values of z and rj_, respectively, at which the expo- 
nent in Eq. (8.20) is w 1/2. Using the expressions for 
C s (r_L,0) and C s (0±,z), given by Eqs. (4.26) (with the 
appropriate wavevector anomalous elastic constants and 
disorder variance) the peak widths (£^) _1 and (£jf) _1 
can be obtained by solving: 



Q- 2 
Q- 2 

which gives 



(«) \ L 



Si 



NL 



Sz S 



G 2 p 



(8.21a) 
(8.21b) 

(8.22a) 
(8.22b) 



The temperature dependence of and could be used 
to determine the exponents r} K , and r)A t since the bulk 
k(T), h(T) in Eqs. (8.21a), (8.21b) have temperature 
dependences that can be extracted from data on bulk 
materials. 
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It is also possible to extract information from the scat- 
tering in the "tails" of the peak. By "tails" one means 
at wavevectors q = Q + <5q such that: 



(£)- 1 «\5q a \«(t; NL ) 



.23) 



where a is either _L or z. In this regime, the scattering 



HSq) = |/>q| 



e -i*q.r e -<(Q.(u(r)-u(0))2>/2 



(8.24) 



is dominated by r < {Sq)' 1 . If \Sq a \ > (£* ) _1 
(the first part of the condition Eq. (8.23)) then the 
exponent ((Q • (u(r) — u(0)) 2 ) is still small, because 
is by definition the length one must go to before 
((Q • (u(r) — u(0)) 2 ) becomes appreciable. This means 
that the exponential can be expanded to give: 



J(5q) 



((Q.(u(r)-u(O)) 2 ) 



1.25) 



which is just the Fourier transform of the quantity be- 
tween angular brackets at wavevector <5q. If |<5<7 Q | <C 
^NLyi ^ thc sccond part f the condition Eq. (8.23)), 

then at this wavevector the fluctuations will be anoma- 
lous and the intensity will be given by: 



[ti(5q)Sql + K(5q)5qj 



(8.26) 



where the l ongitudinal and transverse pieces of 
(ui(q)uj(— q)) have, for simplicity, been set equal, since 
at \6q a \ <C CtvL' tne longitudinal and transverse prop- 
agators are essentially the same (since A(5q)//z(5q) — > 
constant). Thus, by examining the tails of the broadened 
peaks one could directly observe the anomalous elasticity. 



IX. FUTURE THEORETICAL WORK 

There are a number of interesting questions that re- 
main to be investigated experimentally. One of them 
is the problem touched upon in section VIII of apply- 
ing a magnetic field in the opposite direction to the 
spontaneous magnetization to "cancel off" the effects of 
the crystal field. The state thereby produced, though 
only metastable, could, in principle, be exactly, rather 
than approximately, soft. This exact anomalous elastic- 
ity is probably not, however, controlled by the columnar 
elastic glass (CEG) fixed point we have been discussing 
up to now. The reason for this is that the CEG fixed 
point described a model with exact rotation invariance, 
to all orders in the rotation angle 9. In contrast, at the 
metastable fixed point of a model with crystal symme- 
try breaking fields, we have simply tuned H to cancel 
off the leading order in 9 (i.e., the |<9 z u| 2 ) term in the 
Hamiltonian. The full model however, is not rotation in- 
variant to all orders in 9, since, for large 9, the crystal 



symmetry breaking field has a different dependence on 
9 than the H • M term (e.g., cos(0) vs. cos(4(9) for an 
underlying cubic crystal). This lower symmetry of the 
"metastable limit" model allows, e.g., the coefficients of 
the cubic (V • u)|9 z u| 2 and quartic |<9 z u| 4 nonlincaritics 
in the Hamiltonian to be different from one another and 
from those of the quadratic |Vu| 2 terms. This should be 
contrasted with the columnar elastic glass Hamiltonian, 
in which these cubic and quartic terms arise only from 
the square of the nonlinear strain tensor Vij, and must 
therefore have the same coefficients (up to known factors 
of order unity) as each other and the quadratic terms. 

A recent study of this problem 32 has found that the de- 
parture of these cubic and quartic coefficients away from 
equality with the quadratic terms is a relevant perturba- 
tion away from the rotationally invariant fixed point to 
a new non-rotationally invariant fixed point. This new 
fixed point exhibits anomalous elasticity with different 
universal exponents than those given above for the rota- 
tionally invariant problem. Further details will be given 
in a forthcoming publication. 32 

We have focussed in this paper exclusively on static, 
equilibrium properties of (putative) spontaneous vortex 
lattices in ferromagnetic superconductors in presence of 
positional and random tilt disorder. Since, in the pres- 
ence of tilt disorder these lattices are glassy, their dynam- 
ics should exhibit all of the interesting slow phenomenol- 
ogy associated with glasses. The interplay of this slow 
dynamics with the anomalous elasticity we have stud- 
ied here should make the dynamic behavior even more 
fascinating. In particular, the dynamics of depinning of 
these flux lattices, which determines their voltage-current 
(I — V) characteristics, should be quite novel. 

Note added in proof. It has recently been shown 60 
that the elastic theory developed here for spontaneous 
flux lattices in s-wave superconductors also applies, with- 
out modification, to superconductors of arbitrary (e.g., 
p-wave) symmetry. 
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APPENDIX A: LARGE SCALE BEHAVIOR OF 
CORRELATION FUNCTIONS 



Cd ± 



(2tt)^ (ql + l)i 



(A6) 



In this Appendix, we give some technical details on where, in this last equation, the ultraviolet cutoff A = 

the solution of the saddle-point equation (4.20) of the 2i\ja was sent to infinity. Now, equation (A5) can be 

replica formalism of section IV B, which we rewrite here rewritten in the form: 
for definitcness: 



dq z 



TQl 



with [cr](u) = ua(u) — J Q dv a(v). Using the fact that 
dx 3?rV2 



f 



we obtain 



(,x 4 + a 4 ) 2 a 7 



dq z 



a > 



(A2) 



oo [ K qj + Kq\ + [a] (v)] 2 8* 1/4 ( Rq 2 ± + [a] I ' 

(A3) 

and hence, equation (4.20) becomes 
3V2TQla{u) f d d± qj_ 1 



16k 1 / 4 



(Kql + [a](u)) 7/i 



Performing the change of variables 



q± = q " ~K ' 



equation (A4) becomes 

—^-^TQ^Hu)] * =1 

where Cd ± is a numerical constant given by: 
- 1 du (w/«o) 2/e 



[a(u)Y 



Y^t^ T ^(u). (AT) 



= 1, 

Taking the derivative of this last expression with respect 
(Al) to u, and taking into account the fact that [cr]'(u) — 
ua'(u), we obtain: 



7 -2d i s-2d, _ (3s/2cj ± TQI \ 1 



4 w«)i^-li5^^)=- 1 



Solving for [cr](u), we obtain: 

M(«)=(-) f , (A9) 

where we called 9 — d± — 3/2, and where ito is given by: 
3V2c d± TQ 2 



Uq = 



4(7 -2d J _)K 1 / 4 K d ^/ 2 



(A10) 



1. (A4) In terms of [cr](u), the elastic propagator in the presence 



of disorder G(q) is given by: 

1-1 



u 2 ifq 2 + 



' Kql + Kqi\ + J u 2 Kq\ 



(u/u ) 



2/9 



)■ 



+ Kq 4 + {u/u )d 

(All) 



(A5) where, in going from the first to the second line we re- 
placed \<j\(u) by its expression, equation (A9). Now, us- 
ing the change of variables x = u 2 ^ e , we can write: 



dx 



-6/2 



u 2 Kq\ + nq\ + (u/u ) 2 / e 2 J Q x + U /e (Kq 2 ± + K q*) 

= J- \ 2 J?lfl 1 



2' 2' ul/^Kql + nqi)" 



(1.12) 



where 2-F 1 is a hypergeometric function. Replacing this 
last expression back into equation (All), and using the 
following limiting behavior when a — > 0: 

2 F% 1-^,2-^,-1]= r( 2 - -)T( 9 -) a 1 "! + 0(a) 

(1.13) 

we obtain the following long wavelength behavior of the 

dq. 



C(r) 



J (2tt)^ 



correlation function: 



+i 



(1 



G(q) 



6 r(2-|)r(fl/2) 



2-0 w (^i + Kgf) 1+8/2 



(1.14) 



We now arc in a position to calculate the large scale 
behavior of the correlation function G p (r) = ([u(r) — 
u(0)] 2 ), which is now given by: 



)G(q), 
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2T 9 r ( 2-*/2)iW2) fdd± ^ dq * 



1 - e lq r 



u Q 2-9' y ~ ■<-'-*->-' J (2n) d ^+ 1 (Kql + Kqiy+ e / 2 ' 
As we did in the last paragraph, we shall evaluate C p (r±) and C p (z) separately. We have: 



C p (r ± ) 



2T e 



u K 1+e / 2 2-6 
2T 



T(2-e/2)T(6/2) 



TTU K 1+e / 2 



(2n) d J 



dq z 



1 - e iq -' 



2tt (g£ + Kq 2 J K ) 1 + d / 2 



0^A r (!) r (l + f) fd d ^q± 1-e^ 



*4^ VJ 1 2 - 

r(i + f) 



(I) 



(2tt) c 



(3+20)/2 

9± 



(1.15) 



where, in going from the first to the second line we used 
the fact that 

2r ( 5/4)r(f + |) a _ (5+S) _ (ii6) 



rV2 



In(Arx). 



(1.17) 



dx 



2u K 1 / i K 

where, in going from the first to the second line, we used 



(x 4 +a) 1 +f T(l + 6/2) 

For a spontaneous vortex lattice in three dimensions, the fact that r(7/4)r(l/4) = 3ttV2/4. 
d± = 2, 9 = d± — 3/2 = 1/2, and we obtain: 

27T(7/4)r(l/4) f A dq ± 1 - 3 (q ± r ± ) 



C p (r ± ) -- 
C p (z) 



3ttu k 1 / 4 K J q 
d d± qdq 



Jo 2tt 



We now turn our attention to the calculation of C p (z). 
We have: 



2^/^^(l-e-)G (q) , 

rri / f) \ f'°° 

(2-0) r < 2 - e ' 2 ^ L dch (1 " e *" )J ^ 1 dq± W^tlW^ ' 



TTUqK 



(1.18) 



where AT^ = (27r d / 2 /r(<i/2)) is the total solid angle (surface area of a unit sphere) in d dimensions. Now, using the 
fact that 



dqj 



A 



di 



, {ql + Kqi/KY+e/ 2 d ± 

and the limiting behavior when q z — > 0: 



(1+0/2) d± d± 

Ia^J 2jF [ T' 1 + 2' 1 + T' 



A 2 A, 



nq 



4 J ' 



F i A i + * i + ^ _^£i ~ r(i + rf ± /2)r(i + f / Kg 4 ^ 
[ 2 ' 2' 2 ' Kq y~ r(i + e/2) \ka 2 ) ' 



we obtain, in three dimensions (d± = 2): 

/ • c 

Jo 



Ci(x) = ln|x| + 7- 



o ( 9 i + ntf/Ky-w 2 \~) Tz' (L21) 



/* 

Jo 



cos t — 1 



(1.19) 



(1.20) 



(1.23) 



Replacing this last expression into equation (1.18), we 
obtain (here we use the fact that K2 = l/(2n)): 

27T(7/4)r(l/4) l-cos(g^) 
3W io ^ £ ' 

V-Lj^i. (1.22) 



where 7 ~ 0.577 ... is Euler's constant, and the fact that 

A z 

Ci(x) — > when x — > 00, we see that J Q z <ix (1 — 
cosx)/x ~ ln(A z |z|) + 7 for A z \z\ 3> 1, and hence the 
long distance behavior of C p (z) in three dimensions is 
given by: 



Ty/2 



2nu KK 1 / i 



L 



x 



C p (z) 



Ty/2 



where we introduced the cut-off A z ~ l/£ for the inte- 
gration over q z , and where, in going from the first to 
the second line, we used the change of variables x — q z z. 
Now, from the definition of the cosine-integral function: 61 



2-kuqKk 1 / 4 



ln(A|z|) +7 



(1.24) 



as claimed in the text. 
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APPENDIX II: ELASTIC PROPAGATOR IN THE where, in going from the second to the third line, we used 



PRESENCE OF DISORDER 



the following results from projection operators algebra: 



In this Appendix, we show how to derive the expres- 

7ab ( 



sion of the propagator G a a b s (q) of the theory defined by 



the Hamiltonian H 0n of equation (5.12), 

H 0n =J2 I ^(q)<(q)4(-q). 

a.b Jc l 

where 



P Q i 7 (q x )P^(q ± )=P Q i (3 (q ± ), 
Pj 7 (q±)P 7 T £(q±) = P^(q±) , 



(2.8a) 
(2.8b) 
(2.8c) 



( 2 - X ) We therefore obtain for G^(q), in the limit 



n -> 0: 



(2.9) 



(2.2) 



Finding G^(q) amounts to finding the inverse of the 
matrix r^(q) such that: 



G^(q)r^(q) = r^(q)G^(q) = 6 a , S bc . (2.3) 



It is easy to verify that, in the limit n — > 0, G^,(q) is 



G^(q) - Gf (q)i%(<Ll) + G5?(q)P^(qi) , 
where 

Gf(q)=Tr Z 1 (q)5 o6 + A^r Z 2 (q), 
Gf?(q) = TT^q) <5 ab + A t9 2 r T 2 (q) . 

We therefore can write: 

«(q)^(q'))o = (27r) d <5(q + q') G^(q) , (2.11) 



(2.10a) 
(2.10b) 



given by: 



where the subscript in (• • -) indicates that the average 
„ is taken with the statistical weight exp(— f3H nn ) /Zq, with 

G^(q)=r(r- 1 ) Q/3 (q) ( 5 ab + A t(Z 2 [(r- 1 ) ] Q/3 (q). (2.4) Zq = Tr(e X p(-/?tf n)). 



But, since 

r a /j(q) - r L (q)P^(q ± ) + IY^P^) , (2.5) 



APPENDIX III: PERTURBATION THEORY IN 
THE ELASTIC NONLINEARITIES 



we see that 

L/3 



In this Appendix, we fill in a few details on the pertur- 
(r- 1 ) Q/3 (q) = T- L 1 (q)P^ /3 (q 1 _) + r T 1 (q)Pj /3 (q ± ) , (2.6) bativc calculation of Sec. V. We shall start by looking at 

the perturbative corrections to A and /i, before address- 
an d ing the relatively more involved corrections to k and A. 

[( r " 1 ) 2 ] a/3 (q) = ( r - 1 ) a7 (<i)( r - 1 ) 7/3 (q) ; 

= [rZ 1 (q)P a L 7 (q±)+r T 1 (q)P Q T 7 (q ± )] > 
X [r^ 1 (q)P 7 L /3 (q ± ) + r T 1 (q)P 7 ^(q ± )] , 



r^ 2 (q)P^(q x ) + r ? 2 (q)Pj /3 (q ± ) , 



A. Correction to A and /u 

Starting from equation (5.20), and using Wick's the- 
orem, we find that the part of the connected average 



(2.7) (-ffint)o> w hich corrects A and \i is given by: 



(H? nt y 0> [\,n] = Yl I drdr ' [M 2 ^(r)a 7 ^(r')<9 z <(r)a z ^(r)a z <(r')a z ^(r'))^ + 



^ 5 Q ^(r)a 7 <(r') ( d z u a a {v)d z u%{v) (d z u b s (v')Y ) 



,2 \ c 



0> 



+ ^ d a ul{v)d,u\{v') { (d z u%(r)) 2 d z u b (r')d z u b (r>) ) C Q> + 



+ '- 9 Q <(r)a 7 <(r') ( {d z u a p {v)Y {d z u b {v')f ) C Q 



(3.1) 



where we have used the symbol = to emphasize the fact that we only show the part which corrects A and /U. Evaluating 
the above averages using Wick's theorem, we find: 



-7^(H? nt ) c > = \ y" < < < (q)<5T?/(q)«$ < (-q), 



(3.2) 



with the kernel: 
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«t(q±) = ~ff I* >< {yMi GfM) G$(q + q') + Aa5 7 (G^(q') G#(q + q') + G^-(q') G#(q + q')) + 
+ M A<M 7 G$(q') G#(q + q') + A/ifefc G^(q') G#(q + q')} , 

= ~ f f l' 4 X {yM 3 G$(q') G$(q') + M 2 ^ (G^(q') G^ b (q') + G^(q') G^ b (q')) + 

+ fiXqiq-/ G^.(q') G^(q') + X^q a q 3 G^(q') G#(q')} , (3.3) 



where Grg(q) is the propagator of the theory defined by the Hamiltonian Hq h of Eq. (5.12), which can be written in 
terms of transverse and longitudinal parts in the manner: 

G^(q) = Gf (q)P a V<u) + (q)P^(qi) . (3.4) 

where (in the n — > limit) 

Gf (q) = TT^q) 5 o6 + A t q 2 z r- L 2 ( q ) , (3.5a) 

G^ fc (q) = Tr-^q) 6 ab + A t q 2 z r T 2 (q) . (3.5b) 

In going from equation (3.3) to (3.3), we approximated G° b (q + q') ~ G^ b (q'), since we only need terms of order q\ 
in <jr?- (qj_). Now, it is easy to verify, using the algebra of projection operators (Eqs. (2.8a)-(2.8c) of Appendix II), 
that 

G$(q') G$(q') = (Gf (q')) 2 + (d± - l)(Gf (q')) 2 , (3.6a) 

Gg(q') Gf 7 (q') - (Gf (q')) 2 Pf 7 (q±) + (^(q')) 2 ^^) , (3.6b) 

G^(q') G#(q') - (Gf (q')) 2 P 4 L Q (ql) + (G^(q')) 2 P£(q'J • (3.6c) 
The expression of <5r? b (q_i_) therefore becomes: 

«rtf (<u) = 4 1 & x + ?f ) [ ( Gf(q/) ) 2 + (di - 1} ( G ^')) 2 ] + 

+ M V? 7 (G^(q')G^(q') +G^.(q')G^(q'))} , (3.7) 
where we used the following rotational averages of single projection operators (Appendix IV): 



(^(q±)) = |f , 



{P^{<U.)) = ^^5 aP . (3.8b) 

Now, in the second line on the rhs of equation (3.7), using the decomposition (3.4)), we obtain 

G£ 6 7 (q')G#(q') = (Gf (q')) 2 P« 7 (ql)P 4 ■ (ql) + (Gf (q')) 2 ^(q±)^(q±) + 

+ Gf (q')Gf (q') (P^q'jP^q'J + Pj^q'JP^q'J) • (3-9) 

Performing the integration over polar angles, we obtain 



M 2 9^ 7 / ^ 4 (G^ 7 (q')G^(q')+G^(q')G^(q')) 
J a' 



H 2 q±S. 



d ± {d ± + 2) y q , 



" f q' z 4 [2(Gf (q')) 2 + (4 - 2)(G^(q')) 2 + 2d ± Gf (q')Gf (q') 



+ 



- , l ^ 2) jf 9? [4(Gf(q')) 2 + 4(Gf(q')) 2 + 2(^-2)Gf(q')Gf(q') 



Now, using the fact that 
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QiQj = qi p lj(<i±) , 

5«=J^(qi)+if(qi) 



we obtain that 5Tf b (c[±) can be written in the form: 



with: 



1 /A 2 2A/z\ 



(Gf(q')) 2 + (^-l)(Gf(q'))- 



+ 



Td ± (dj 



-L J q' J 

^q^y [6(Gf(q')) 2 + (24-2)(G^(q')) 2 + 4(d ± -l)Gf(q')G^(q 







srf ( q± ) = - 



Td ± (d ± + 2) 



/q' 



2(Gf (q')) 2 + (4 - 2)(G^ b (q')) 2 + 2d ± G? (q')Gf (q' 



(3.10) 



(3.11a) 
(3.11b) 



Now, <5r^ b (q^) is nothing but the correction to the 
transverse part of the elastic tensor Bt, i.e. ST^'(c[±) = 
<5/i. On the other hand, 5T'^ b (q±) is the correction to the 
longitudinal part Bl = A + 2^. Hence: 

(5A = <5(A + 2/i) - 25^ 

= fff (qj.) - 25lf (q ± ) (3.12) 

Performing the above q- integrations, we find that both 
corrections S/j and SX behave like 

5n ~ 5\~ ^^K^At L 7 - 2d (3.13) 

where L z is the size of the system along the direction 
of the flux lines (which comes in as a result of imposing 
an infrared cut-off 1/L Z on q z integrations), which is the 
result (5.22) quoted in the text. 



5D = gDF D (x), 
Sfi = g/j,F^(x), 



(3.14a) 
(3.14b) 



where x = X/fi, and where the functions Fd and F^ are 
the functions defined in Eqs. (6.7a)-(6.7b) of the text. 



B. Corrections to k and A t 



We now turn our attention to the calculation of the 
perturbative corrections to the curvature modulus k and 



In the RG analysis of Sec. VI, the integrals in Eqs. to the disorder strength A t . These corrections will 
(3.11a)-(3.11b) are evaluated within the momentum shell come from the following terms of the connected average 
defined by Eqs. (6.1a)-(6.1b), and this leads to: (Hint,n)o>> 

(C,n)o>[«-At]=E / drdr'{ M 2 a zU ^<(r)a zU b <(r')(a Q ^>(r)a z ^>(r)9 7U b > (r')d z u b > (r')) 

a, b 

+ tfd z u a < (r)d z u b < {v'){d a uf (r)d z uf> (r)d 7 u b> (r')d z u b > (r')) c 



o> 



'o> 



+ tfd z u a < (r)d z u b < {v')(d a uf (r)d z u a > (r)d,u b > (r')d z u b > (r')) c Q> 
+ tfd z uf (r)d z u b < (r'){d a uf> (r)d z u a > (r)d,u b > {v')d z u b > (r'))^ 
+ n\d z u a < (r)d z u b < (r')(d a u a > {v)d z uf (r)d,u b > (v')d z uf (r')) c Q> 



+ ii\d z u a < {v)d z u b < {v'){d a uf (r)d z u a > (r)d 7 u b > (v')d z u b > (r')) 



o> 



X 2 d z up {v)d z u b <{v')(d a u a >{v)d z uf (r)9 7 <>(r')0^ b >(r'))^}. 



(3.15) 



All the connected averages on the right hand side of the above equation are of the general form 
(9 Q u^ > (r)9 z Mp > (r)9 7 'u^ > (r')9 z 'u b> (r')) 0> . Using Wick's theorem, we can write: 



{d a uf (r)d z u a > (r)9 7 u b > (r')d z u b > (r')) c Q> = 



q^q'lGfMGfM) + q a q z q' z q\GfM)G a p h S (ci') 



*(q+q')-( r - r ') 



Integration over r and r' leads to the result 
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J drdr^ U J<(r)0 z i£<(rO(0 Q u;^^ x 

x ^ \q a q' 7 (q z + <£) 2 G$(q')G£(-q - q') + q a q' z (q z + <£)(g 7 + g^G&faOG^-q - q') 

Jq' L 



(3.16) 



Using the above general result, we calculate the seven 
terms appearing on the right hand side of equation (3.15), 
and then Taylor expand the result in q z near q z — 0. Te- 
dious but straightforward calculations (which are best 
carried out using a symbolic math processor such as 
Mathematica) lead to the result: 

~(Hf nt )o> = |%r(q)«t(q) Uj fc< (-q), (3.17) 



where 8Tfj(q z ) is given by 



Sr%(q z ) = S K q% b - 5(^)ql8 %] . (3.18) 

Again, it can be shown that evaluating the q' integrals 
in Eq. (3.16) with an IR cut-off 1/L Z on q z integrations 
leads to same conclusion as in the previous paragraph, 
namely that 5k and 5(A t /T) both diverge with the sys- 
tem size in the fashion 



5k ~ S(A t /T) ~ L 



7-2d 



(3.19) 



On the other hand, if the q' integrals on the rhs of Eq. 
(3.16) are evaluated within the momentum shell of Eqs. 
(6.1a)-(6.1b), then it follows that: 



5k = KgF K (x) , 
5^)=A t gF At (x). 



(3.20a) 
(3.20b) 



where the functions F K and F At are the functions defined 
in Eqs. (6.7c)-(6.7d) of the text. 



APPENDIX IV: ROTATIONAL AVERAGES OF 
PROJECTION OPERATORS 

In this Appendix, for completeness we evaluate angu- 
lar averages of projection operators, necessary for the 
RG computations in the main text and Appendix III. 
More explicitly, we want to express the following d±- 
dimcnsional integrals: 



1 a/3 



rLL 



[/]= / /(3i)i^(<u), 

Jq± 

[/]= / /(«i)i^(<u), 

Jq ± 

[/]= / f{ql)P^) P , L M^ 

Jq i 



(4.1a) 
(4.1b) 
(4.1c) 



Ills\f]= ! /( < zl)Pj /3 (q i )P 7 i 5 (q i ), (4.1d) 

i«e lS [f]= I f(ql)PU^) p ^\ (4.ic) 

Jq± 

where f{q\) is an arbitrary function of q\, in terms of 
the spherically symmetric integral 



I[f] 



f{ql). 



(4.2) 



We start with the integrals I%p[f] an d I^p[f] which can 
easily be evaluated. We have (in what follows, q a denotes 
the quantity q a /q)' 

Jq± 

= 5 afj [ f(ql) q l, (4.3) 

jq ± q 

where no summation is implied on the index a in the 
second line. Using the rotational symmetry of /, we can 
write: 



l a/3 



i.e. 



(4.4) 



(4.5) 



Similarly, we have for Jj^: 

Q[f] = I /(«i)[«a/J-i^(qi)], 

Jq i 



Mi4' / /( ^ } < (46) 

d± Jq ± 



so 



(4.7) 



We now turn our attention to the integral I^p -ys[f]- 
We first consider the case where the indices are all equal, 
i.e. we want to calculate the integral: 
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J a i 



L 



qi 



7T fill) ■ 

qi Ij- 



(4.8) 



(ji = q sin9d ± -i sin 9 d± - 2 1 1 1 sin 02 cos 6*i 
q 2 = q sin _ i sin _ 2 • • • sin 9 2 sin #i 
q 3 = q sin9 d± -i sm9 d± - 2 • • -sin^ 3 cos# 2 
<74 = q sin 0^ - 1 sin 9 d± - 2 • • • sin #4 cos 9 3 

Qd±-i = 1 sin6» di _i cos# di _ 2 



and hence we see that the most convenient choice of a for 
calculating the integral (4.8) is a = d±. Taking q a = q d± 
results in the following expression 



Due to the rotational symmetry of the function /, this 
integral will be the same for all possible values of the in- 
dex a = 1, • • • , d± . The decomposition of the vector qj_ 
in spherical coordinates in d± dimensions is given by 



r" [/]= f^SL 



d d± q± = q d± - 1 sm d ^- 2 9 d± - 1 de d± - 1 sin^" 3 
(where < 9k < tt for k ^ 1, and < 9\ < 2n), we obtain: 



(2TT) d 

Using the fact that 

• sin #2 d9i d9\ 



cos 4 6^_i f{q\). 



(4.9) 



(4.10) 



2 P PIT PIT />2tT 

I«a, a *[f] = j^yTZ J <Zx ±_1 /(<7x) dq± J sin^- 2 6 d± -i cos 4 e d± -! d9 d± . 1 ...j sm9 2 d9 2 j d9 1 , 
= C d± I f{ql) 

J a 1 



(4.11) 



where 



a 



d± 



Jo Sm 



i 1 -2 



9 cos 4 9 d9 



j* sin^" 2 6 d9 



(4.12) 



The integrals in the numerator and denominator of C d± 
can be evaluated analytically. Mathematica for exam- 
ple gives: 



C d± = 



r(d±/2) 



30? 



V r(i 



r(d±/2) 

4r(^) 2r(^L) 
r(2 + d ± ) + r(3 + d±) 



3r(^) 

(l + dx) 

) 



. (4.13) 



Expanding this last expression using the following results 
for the Gamma function: 



T(l + z) = zT(z), 



(4.14a) 
(4.14b) 



leads, after a few manipulations, to the very simple re- 
sult: 



C d± 



d ± (d ± +2) ' 



Hence: 



tLL 



d±(d ± +2) 



I[f\- 



(4.15) 



(4.16) 



We now consider the case where the indices a, (5, 7 
and 5 in Iap ~g are not all equal. In that case, we can 
write: 



rLL 

-'a/3,75 



= / /(9i)^(<u)^ 4 (cii), 

Ja± 



f(ql)q a q(3q~(q5 , 



-/ 

= S a pS 7S / f(q 2 ± )qlq^ + S ai Sp S / f(qi)qlqj 

+ S aS 5 Pl [ f{ql)qlq} , (4.17) 
J qj_ 

where the two remaining indices in each integral are dis- 
tinct from each other. Let us evaluate one such integral. 
We have: 

/ f(q 2 )q 2 jj = [ 



2 2 

f(q±)—r- > 



dx- 



dx 



dx 



1 Jq ± q± 

i —[ f if&)-[ 

1 Jq, 9 I ./q, <?X 



(4.18) 



But: 
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/ f(il)$= I /(<d)# = ^-/ /(£). (4. 

./q , (/_L V «-L Jq , 



'q_L ?_L -/q_ 

and 

„4 



19) 



Jqj_ yj_ Jq^ Jq^ 



so we obtain: 



Using the result (4.15) for Cd ± , we finally obtain: 

Ij^-u^mlj^ ,4 - 22) 

and hence, Eq. (4.17) gives (a, (3, 7 and 5 not all equal): 

(SafjSjS + SayS/SS + Sasfip-,) I[f] ■ 

(4.23) 



d±(d± + 2) 



It is not difficult to see that equation (4.16) is a spe- 
cial case of this last equation. Hence, the final result 
for Iaf> tl s[f] is j us t equation (4.23), which describes all 
possible combinations for indices. 

Now, let us find I^slfl We have: 
I*l~ fS [f}= I /(«±)^(<ll)££(<U), 

/(gi)P^(<Ll)[5 7 «-P>U)], 

= j-^Mf] d±{d ]_ + 2) x 

x (5 a ()8 1 s + 5 ai 5ps + SasSfjy) j /[/] , (4.24) 



hence 



- d±id l + 2 ) ( 6a ^ 5 + 5asS f>J } I[f] ■ (4 - 25) 

The last rotational average that we need to calculate is 
lal-yslfl We have: 

iZLslf}= [ /(«i)i^(<u)i^(qi), 

= / /(<?!)[ <W - i^(qi)] [<M - P 7 L 5 (q±)] , 
q 1 



= <5 7 <5 / /(<7i) [8 a p8 7S - <5 Q/3 P^(q_L) 

- 5 7 «i^(qi) + i%(<ll)i$(<Ll)] , 
2 

= ^a/3^7«5 - 3— ^a/3^ 7 5 -f [/] 

+ d ± (dx + 2)^ 4 + 
+ S ai Sps + S a sSp 7 ) I[f] , 



(4.26) 



hence 



(d± + 2) 



^a/3<5 7 (5 + 



+ 



(<W/3* + <W^ 7 ) } /[/] • (4-27) 



d±(d ± + 2) 
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